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PREFACE 


This issue of the Naval Research Logistics Quarterly contains the 
majority of papers which were presented at a Conference on Com- 
putational Aspects of Linear Programming sponsored by The RAND 
Corporation on 30-31 August 1956. The purpose of this Conference was 
to ''present a broad survey of the status of computation techniques for 
those who formulate models or administer code developments.'' The 
continual development of more powerful computational methods for the 
solution of linear programming problems enlarges the scope of its pos- 
sible applications. While the mathematical disciplines involved are of 
interest primarily to the specialist in the field, significant advances in 
this as well as other mathematicalareas have an important meaning for 
the nonspecialist working toward a more rational basis for logistics 
systems. It is hoped that the publication of the proceedings of the 
Conference in this issue of the Naval Research Logistics Quarterly will 
provide a useful function in disseminating these papers to a wider 
audience of specialists and nonspecialists alike. 


The organization of the papers in this issue does not follow the 
order in which the papers were presented at the Conference. In the 
interest of accurately recording the Conference these papers have not 
been subject to editorial review, contrary to the established practice of 
this journal. 


Abstracts for the following papers were not received in time for 
inclusion: "Types of Computing Jobs and Red Tape" by W. Orchard - 
Harp (RAND); "Integrated Logistic Management of an Air Force Item" 
by J. R. Gillespie (Air Materiel Command, U.S. Air Force); and "Salvaging 
Techniques after Data Error" by L. Cutler (RAND). 


The editors wish to express their appreciation to L. R. Ford, Jr., 
D.R. Fulkerson, and G. B. Dantzig for their cooperation inthe publication 
of this issue. 


—M. E. R. 








PROCUREMENT AND ALLOCATION OF NAVAL 
ELECTRONIC EQUIPMENTS! 


George Suzuki 


Bureau of Ships 


In its simplest conceptual form, the allocation problem I wish to describe is this: 

Each naval ship has certain functions which it must perform with electronic equipments. 
We shall call a specific combination of a function and a ship an "activity." Each activity, j, is 
authorized to receive up to its allowance of equipment sets, b., and although a set may be 
capable of performing more than one function, a set once installed is assumed to fill only one 
allowance position. There are m equipment models and each model, i, is available for assign- 
ment in quantities a;- The problem then is, how should these available sets be assigned to the 
fleet in the best way ? 

Suppose we hypothesize a set of values, Cijs for each possible assignment situation such 
that these values represent in an acceptable manner the military worth of the assignment of an 
equipment model i to activity j. These values are determined such that Ci > Cen implies that 
assignment of one set of model i to activity j is preferable to the assignment of one set of 
model k to activity h and furthermore, such that Ci, = WCpp, W > 0, implies that the assign- 
ment i,j is w times more desirable than assignment k,h. With such values, we can formulate 
the allocation problem as a linear programming problem of the transportation type and deter - 
mine allocations which are optimal within this framework. That is, determine Xj; 2 0, 
i=1,2,...m,j=1,2,...n, where Xj represents the quantity of model i allocated to activity 
j, such that 


: ra Xij is maximum 


subject to the availability restraints 


ct b 
z Xj = j 





l cos - , 

For general description of this problem see Jack W. Smith, ''A Plan to Allocate and Pro- 
cure Electronic Sets Using Linear Programming Techniques and Analytical Methods of Assign- 
ing Values to Qualitative Factors,'' Naval Research Logistics Quarterly, Sept. 1956 
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The procurement problem is, given a sum of money, M, and unit prices for each model, 
Pi, how should the money be allocated to the equipment models to effect an optimum purchase? 
A good procurement plan must, of course, consider equipment in stock or on order and the 
manner in which the purchased equipment is to be used by the fleet. Thus procurement and 
allocation must be considered simultaneously. We have formulated the problem as a linear 
program as follows: 


Maximize ; 5H Xj 


subject to the restraints 


Toa. 5%, ’ 
Ey 5% * % 


The y;'s represent quantities of each model to be purchased. 

In this more general formulation, if we insisted on obtaining an optimum integral valued 
solution, the computational problem would be quite difficult. At least Iam not aware of any 
straightforward systematic method of obtaining such a solution. However, the restrictions of 
this problem are such that when solved by the Simplex Method, there can be only a few non- 
integral values in the optimum solution. Because of this, we have used the simplex method in 
solving the procurement allocation problem and have selected the nearest integral coordinate 
point by rounding the nonintegral values. 

The individual problems we deal with in the Bureau of Ships do not include all types of 
electronic equipments and all possible activities simultaneously. From the procurement point 
of view, this general problem is more useful. However, the difficulties encountered in such an 
overall formulation have not been resolved. The problems of the type I wish to describe are 
subproblems involving individual categories of equipment applications. It is assumed that the 
total equipment budget has been allocated to each equipment application category and that these 
sums are given. Each category of equipment application is then treated as an independent 
problem. To illustrate, the applications of electronic equipments are stated in categories like 
air-search radar, or radio transmitter, or echo sounder, or surface-search radar, or radar 
repeater. Each category is divided into subcategories or requirements. For radar repeaters, 
the requirements are standard watertight repeater, standard nonwatertight repeater, and off- 
center repeater. For radio transmitter, the requirements are 0.3 to 2 mc, cw/voice, 50 miles 
range; 2 to 26 mc, cw/voice, 200 miles, etc. A single problem consists of only those require- 


ments within a category which have common equipment model applicability and those equipment 
models which have requirements in common. 
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The usefulness of the linear programming solution to this problem obviously depends 
upon the values assigned to the C4;'8 and so I shall briefly describe their derivation.2 It must 
be emphasized that the values we obtain are not an attempt to determine, or even approximate, 
real worth since there is, of course, no absolute numerical measure of "military worth." The 
values we derive, however, are numerical values which are consistent with certain military 
opinions and decisions with respect to desirability of replacements and importance of activ- 
ities. The procurement-allocation plans we obtain are at least consistent with these opinions 
and decisions. 

Each activity generally has equipment already installed so the new equipment which is 
assigned will be a replacement. It is necessary to consider the installed equipments as a 
factor in determining assignment preferences since, other things being equal, it is clearly 
more desirable to replace a poor set than a good one. Each activity thus is considered to have 
an associated installed equipment (or equipments in case of multiple allowances). In case an 
allowance is not filled and there is no installed equipment, this state is called a void, that is, 
a "void" is equivalent to a poor equipment model. Each unique activity-installed equipment 
combination will be called an "activity state.'' We now change our index designation slightly 
from the initial discussion so that the c;,'s represent the military worth of assigning one set 
of model i to activity state j, where activity state j represents a unique combination of ship 
type, requirement, and an installed equipment model. 

Each Ci is considered to be a function of two factors. (1) The improvement to be 
gained by assigning an equipment model A with "goodness" value A., to replace an installed 
model B, with "goodness" value Bj, in the activity state j; and (2) the military importance of 
the activity (ship type and requirement) corresponding to activity state j. We hypothesize that 


tg? Ay SPs 


where the index i corresponds to model A, (A; - B) represents the improvement and m, the 


J 
importance of the activity. 


The goodness (or suitability, or preference) of an equipment, which includes the evalu- 
ation of maintenance difficulties as well as performance, depends upon the activity to which it 
is assigned. That is, equipments capable of fulfilling several requirements may not fulfill each 
requirement equally well and equipment size or design may not be compatible with the struc- 
ture of certain ships. Thus, the goodness values of the equipments must be established for 
each activity. In practice, we have set up activity classes composed of several ship types. 

This has considerably reduced the amount of work to be done, but for this discussion, the 
description will be in terms of individual ship types. To determine the goodness values by our 
scheme, we ask authorized naval officers questions which will yield the following type of infor- 
mation for each activity: 

(1) Equipment models in order of preference. Say, A>B>C>...>N. 

(2) Positive differences between models in order of preference. For example: 

A-G>B-G>A-F> etc. 
This is interpreted as: 





2A complete description of the derivation by Robert J. Aumann and Joseph B. Kruskal will 
appear in a future issue of Naval Research Logistics Quarterly 
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Model A replacing Model G is preferred over Model B replacing Model G, 


Model B replacing Model G is preferred over Model A replacing Model F, 
etc. 


From this information we derive numerical values for each goodness which will be consistent 
with the preferences indicated. To illustrate, suppose there are four equipment models 
ordered as 


(1) A>B>C > Void 
for a particular activity. Suppose the differences are ordered as 


(2) A - Void > 
B - Void > 
A-C> 
C - Void > 
B-C> 


To determine numerical values, let 


(3) 


Then by (3), (1) and (2) are equivalent to 


(4) a+b+c=l 
a+b>c 

b>a 

e>p 


a>0O 


Generally, there is no unique solution to the equalities and inequalities resulting from the 
ordered differences. At present we have no basis for preferring any particular solution. In 
the absence of any such criterion, we have used the centroid of the polyhedron described by the 
system, as illustrated by (4). To determine the centroid, we determine all of the vertices of 
the polyhedron and take the average of the coordinates of the vertices. Since there are rela- 


tively few equipment models in any problem, the determination of all vertices is not difficult. 
For (4) the three vertices are 
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a=1/3 a=0 a=1/4 
at tas b= 1/3 V9 = b= 1/2 ~—* b= 1/4 
c= 1/3 c=1/2 c=1/2 


and taking the linear combination 


v,/3 + v3/3 


the centroid is 


This yields a set of goodness values 


A=1 
B= 29/36 
C = 16/36 
D=0 


The smaller the polyhedron, the more desirable, since this would reduce the possible 
variations in the procurement allocation plans. At times, the decisions will yield equalities 
rather than inequalities. Since each equality reduces the polyhedron by one dimension, naval 
authorities are urged to indicate equalities rather than to indicate inequalities based on fine 
hair-splitting decisions. 

The goodness values are determined relative to the most preferred equipment model 
for each activity. Since equipment models may be applicable to several requirements, hence 
activities, with varying goodnesses, the goodness values for the activities should be related. 
That is, the goodness values of equipment models applicable to several activities must be 
consistent with the preferences between activities as well as with the other models applicable 
within each activity. The procedure we have used so far for doing this is similar to that of 
obtaining the importance values, m,, which will now be described. 

The activities are ranked in priority order by the Chief of Naval Operations (CNO) in 
the form of a Military Improvement Plan (MIP), so that we are given, say, 


m,; ~My > Mg >... M,.- 


The rankings are established independently of equipments associated with the activities and 
constitute CNO's opinion of the importance of the activities. To establish the values for the 
m's, we assign 


m,=1 
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and the remaining values are determined relative to m1; by a series of questions of the fol- 
lowing sort. 

"Suppose one equipment model A is available for assignment and there are two activity 
states to which the set can be assigned, 


activity state j (with importance m,) with equipment model B installed 
and activity state k (with importance m,,) with equipment model C installed. 
Where would you assign the available set A?" 
The answer will imply that 


(A, - B,) m, > (A, - C,) m, 


(A; - B;) m, < (A, - C,) m, 


Since numerical values for (A; - B;) and (A, - C,) are available, the preference yields 


ms A; -B 


eo Se | 


m, A, - By 


7. 


m, Ay - Cy 


By having the goodness difference ratios in numerical sort and by changing the equipment 
models in the questions, it is possible to ask questions until an equality is attained or a change 
in the direction of the inequality results. This establishes the importance ratio to a particular 
value or to an interval. In practice, the difference ratios appear to be sufficiently dense so 
that the intervals established for the importance ratios are small. 

In practice, we have not determined the importance values for each activity. There are 
priority classes within the MIP structure where major differences exist between classes and 
differences within a class are relatively small. We determine importance values for classes 
of activities and assume that the importance values are the same within a priority class. The 
priority class is, like the activity class, a group of ship types associated with specific require- 
ments. The ship types for the priority and activity classes need not be the same since the 
criteria for classification are different. The criterion for the priority classes is primarily 
the mission of the ships while for the activity classes, it is primarily the structure of the 
ships in the sense of feasibility of installation of equipments. 

The scheme as described has been tried on one problem involving surface-search 
radar. The results appear to be satisfactory in the sense that no improvement in the 
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procurement -allocation plan yielded by the scheme has been forthcoming. Many problems 
remain to be resolved before the scheme can be considered for routine use. However, the 
procedure seems to offer possibilities of injecting objectivity in an area heretofore dominated 
by intuitive and qualitative methods. 





LINEAR PROGRAM PLANNING OF REFINERY OPERATIONS 


L. K. Cheney 
Richfield Oil Corporation 


Many complex technical economic problems encountered throughout the entire petro- 
leum industry may be effectively evaluated by use of linear programming techniques. This 
method, especially when coupled with the use of modern high-speed digital computers, allows 
the solution of very complex problems which hitherto were unsolvable. 

Typical of these problems is a relatively simple gasoline-blending, catalytic -reforming 
and finished-product-mix model which Richfield has employed to guide its refinery operations. 
The matrix of this problem contained 27 equations and 66 variables. 


THE PROBLEM 
Given: 

1. Thirteen gasoline-blending components and one kerosene-type stock along with 
quantities and properties of each. 

2. One catalytic-reforming unit with upper and lower capacity restraints, catalytic- 
reforming yield curves for varying severity of operation, and catalytic-reforming cost curves 
for varying severity of operation. 

3. Excess normal butane available in unlimited quantities, at a fixed cost, necessary 
to blend the finished products to a specified vapor pressure. 

4. Unlimited tetraethyl lead as required, at a fixed cost, to blend the finished products 
to specified octanes, but with an upper limit on the use of tetraethyl lead in gasoline blending. 

5. Aclose approximation of the effect of tetraethyl lead in the finished blend of each 
product with respect to octane generation. 

6. A fixed quantity of premium and regular gasoline to be produced, no restraint on 
the quantity of jet fuel (JP #4) to be produced, an upper limit on the amount of military combat 
gasoline to be produced and no restraint on the amount of kerosene-type stock produced as a 
product. 

7. Two octanes, one volatility and a fixed vapor-pressure specification for premium 
gasoline; one octane and a fixed vapor-pressure specification for regular gasoline; two octanes, 
one volatility, and a fixed vapor-pressure specification for military combat gasoline; and four 
separate specifications for jet fuel (JP #4) product. 

8. Sales realization value for each finished product. 
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The optimal premium gasoline blend, regular gasoline blend, military combat gasoline 
blend and quantity, jet fuel (JP #4) blend and quantity, kerosene-type stock product quantity, 
and catalytic-reformer charge rate and severity of operation to obtain the maximum reali- 
zation from sales of the products with proper debits for costs incurred. 

Figure 1 is a simplified block flow diagram of the problem. Note that the input con- 
sists of 14 base stocks plus tetraethyl lead. One of these base stocks, straight-run gasoline, 
may be used directly as an input or it may be further processed through a catalytic -refor ming 
unit to produce reformate which then becomes a new input. Output consists of four finished 
products plus a portion of one of the unused input stocks, kerosene distillate, plus the by- 
product produced from the catalytic reformer. 
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Figure 1 


























Typical quality specifications for one of the output products and typical quality inspec - 
tions for three of the input base stocks are given in Table 1. These data are representative of 
similar information required for all input and output streams. 


TABLE 1 
Typical Inspection and Specification Data Required for Input and Output Stocks 





Research | Motor | %Evap. | Reid Vapor 
Octane /| Octane | at 235°F | Pressure 








Premium Gasoline Specification 97.0 86.0 50.0 9.5 
STOCKs* 


NC, - NC, 78.6 74.9 100.0 9.5 


St. Run Gasoline 80.9 76.5 33.1 9.5 
90 RON Reformate 96.4 87.4 43.4 9.5 























*Stocks are preblended with butane to 9.5 RVPand are ethylized to 2.0cc/gal TEL. 
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Figure 2 represents the effect of varying 
quantities of tetraethyl lead on the research 
octane number of the final premium-gasoline 
product whichis anticipated. Similar curves are 
required for each type of octane for each kind of 
product. The validity of the final linear program —FESEARCH 
solution is dependent on the reliability with which a 
these various functions can be predicted. This 
type of function is obviously nonlinear so that it 
must be built into the linear programming model 
by use of linear approximating procedures with 
proper upper bounding equations. In practice this 
particular functional representation does not 
require lower bounding equations because of its ‘ 
convex nature. The use of tetraethyl lead to TETRAETHYL LEAD, CC/GAL 
produce octanes is costly and since the first 








Figure 2 - Premium-gasoline octane 


cc/galof leadgives agreater increase in octanes vs TEL 

than any subsequent cc/gal, the linear program- 

ming iterative procedure will automatically use 

the lowest lead content possible. The error introduced by the linear approximation procedure 
can be decreased to as small a value as desired by picking increasing numbers of lead incre- 


ments. In the example given in Figure 1 it can be observed that if the final solution falls 
within the 1 to 3 cc/gal lead level, the approximation is well within experimental error of 
octane measurement. However, if the solution should fall between 0 and 1 cc/gal of lead, 
smaller increments should be designated. 

Figures 3 and 4 represent basic-yield and cost information as a function of the severity 
of catalytic-reformer operation as measured by the clear octane of the reformate produced. 
These functions are similar to the lead-octane curve of Figure 2 in that they are nonlinear. 
Fortunately these functions are also convex with respect to yield and cost relationships so 
that the final solution of the model will automatically pick the lowest severity consistent with 
the restraints of the problem. The higher the severity of reforming required, the more costly 


REFORMATE YIELD FUEL GAS YIELD 
VOLUME % 90 ate VOLUME % FOE 











90 
REFORMATE RESEARCH OCTANE , NEAT 


Figure 3 - Catalytic-reforming yields vs reformate octane 
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it will be both with respect to yield of 
products produced and the cost of oper- 
ation. These functions were also built into 
the model by the linear approximation 
INCREMENTAL technique. However, as will become 

~—_— apparent, they were handled by multiple- 
activity introduction rather than by 
increasing the number of equations in the 
model as was done with the lead-octane 
relationship. 


cosT 
$/88L FEED 








REFORMATE RESEARCH OCTANE , NEAT THE MATRIX 
Figure 4 - Caialytic-reforming costs vs The major effort required in the 
reformate octane solution of any complex problem by the 
linear programming procedure is in the 
formulation of the problem into a meaning- 
ful matrix. The value of the solution obtained is completely dependent on the ability to reduce 
the relationships to the matrix form in a manner consistent with the physical facts. A repre- 
sentative portion of the above stated problem is presented in Table 2. The activities and 


equations given were selected to illustrate the diverse nature of both the equations and activi- 
ties which were employed. 


TABLE 2 
Representative Portion of the Total Problem Matrix 





Activities 1 2 3 4 5 


NC;-NCe , 90 RON 
Equations To To Reformate 
Prem To Prem 





Objective 5.63 ¥ 5.08 





NC,-NC¢ Availability 1.026 
St. Run Availability 
Premium Quantity 
Premium % at 235 














Premium RM Octane 





Premium MM Octane 





Premium TEL Balance 



































Activity (1) 

This is the activity of routing the amount of NC -* NC, base stock pressured to 9.5 
Reid Vapor Pressure with excess butane and leaded to 2 cc/gal to produce one barrel of 
premium gasoline. This activity uses 1.026 barrel of NC; - NC¢ stock availability, produces 
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one barrel of premium gasoline, exceeds premium-gasoline volatility specification by 50 units, 
is short premium-gasoline RM octane specifications by 18.4 units, is short of premium- 
gasoline MM octane specifications by 11.1 units and makes available one barrel of premium 
gasoline which can be further leaded for octane appreciation. It also makes available, for 
subsequent use, 0.026 barrels of butane since the vapor pressure of the NC .* NC, stock 
actually exceeded the premium-gasoline Reid Vapor Pressure specification of 9.5. 


Activity (2) 

This is similar to activity (1). It uses 0.847 barrels of straight-run gasoline stock 
availability, produces one barrel of premium gasoline, is short of premium-gasoline volatility 
specification by 16.9 units, is short of premium-gasoline RM octane and MM octane specifi- 
cations by 16.1 and 9.5 units respectively and makes available one barrel of premium gasoline 
which can be further leaded for octane appreciation. It also utilizes 0.153 barrels of excess 
butane to bring it to the premium-gasoline vapor-pressure specification. 

In the total problem, activities such as (1) and (2) were introduced for all input stocks 
except kerosene distillate to produce all products except jet fuel (JP #4). 


Activity (3) 

This is the activity of processing straight-run gasoline stock through the catalytic 
reformer to produce net by-product and 90-octane unleaded reformate which is pressured to 
9.5 Reid Vapor Pressure and leaded to 2 cc/gal to produce one barrel of premim gasoline. 


This activity uses 0.977 barrels of straight-run gasoline stock, produces one barrel of pre- 
mium gasoline, is short premium-gasoline volatility specifications by 6.6 units, is short of 
premium-gasoline RM octane specifications by 0.6 units, exceeds premium-gasoline MM 
octane specifications by 1.4 units and makes available one barrel of premium gasoline which 
can be further leaded for octane appreciation. 

In the total problem, activities such as (3), at 85, 90 and 95 unleaded octane severity, 
were introduced to produce all products except jet fuel (JP #4). This is an example of using 
multiple activities to linearly approximate a nonlinear function, as previously discussed. The 
optimal solution could result with one or more of these activities in the basis. This would 
mean the solution would require an operation between these activities, the exact point depend- 
ent on the ratio of one to the other. 


Activity (4) 

This is the activity of adding lead of the type associated with increasing lead content 
from 2 to 3 cc/gal in premium gasoline. This activity increases the total premium-gasoline 
RM octane by 0.8 units, increases the total premium-gasoline MM octane by 0.9 units, and 
utilizes one barrel of premium gasoline which is available for further lead additions. 

In the total problem, activities such as (4), for removal of lead of the type associated 
from 2 to 1 cc/gal and from 1 to 0 cc/gal, are introduced for all products except jet fuel. 
Associated with each activity of this type is one additional equation limiting the amount of 
product which is available for each type of lead change. Again this is an example of linear 
approximation used to approximate nonlinear functions. As previously discussed this proce- 
dure requires not only an addition to the number of activities, but also requires the addition 
of upper bounding equations. In this case, the optimal solution could result with a basic 
Solution leaving a positive value in the slack of one of the upper bounding equations. The actual 


427791 O- 57-2 
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amount of lead to be utilized can be determined by the ratio of slack to the total available 
quantity of the product which could have been used for lead addition of that same type. 


Activity (5) 

This is the activity of producing one barrel of jet fuel (JP #4). It utilized 0.103 barrels 
of NC .* NC, stock availability and 0.897 barrels of straight-run gasoline availability to 
produce one barrel of specification jet fuel. 

In the total problem six additional activities representative of several jet fuel recipes 
were introduced. These recipes varied greatly in composition. They were designed to be 
limited by one or more specifications, each activity representing a different limitation. The 
intent was to include for all practical purposes nearly all the possible combinations within the 
theoretical convex. This procedure is dependent on the problem formulator's ability to closely 
approximate the overall convex, but it is particularly useful where there are many stocks but 
a relatively few specifications since it adds activities but not equations to the matrix. The 
optimal solution may contain several of these activities in the basis, but a combination of 
these activities represents the final product blend and the meeting of specifications is assured. 


SLACK 

For purposes of presentation the slack variables of the partial matrix given in Table 2 
are listed in one column. It should be understood that the slack variable of each equation is 
a different variable. All of the equations except the premium-gasoline quantity and the 
premium-gasoline RM octane have associajed slack variables. These equations have no slack 
since the problem requires these quantities to be met exactly. A similar criteria was observed 
for all equations in the total problem. 


OBJECTIVE 

The objective coefficient for each activity may be defined similarly for all activities. 
The problem was set up in the form of a maximized objective. 

The objective coefficient for each activity equals the sales realization from each 
product produced, less the value of butane utilized, less the cost of the tetraethyl lead pre- 
blended with the component, less the operating expense incurred if further processing was 
done, plus the value of by-products produced by any required further processing. The one 
exception to this basic definition was for activities like Number 4, where the objective form 
is defined as minus the cost of adding the amount of lead specified to one barrel of product. 


RIGHT-HAND SIDE 

The right-hand side of the matrix represents the quantities to which each of the equa- 
tions are equal. Thus the availability of NC “ NC, stock was 3000 barrels. The availability 
of straight-run gasoline stock was 14,000 barrels. The quantity of premium gasoline required 


was 25,000 barrels. The sum of the negative and positive contributions to the premium- 
gasoline volatility specification shall be zero etc. 


THE SOLUTION 

The solution to this particular problem was obtained by using the IBM 704 digital com- 
puter. From the optimal solution, the availability of military combat gasoline, jet fuel, and 
kerosene distillate available for sales was determined after setting aside the desired quantities 
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of premium and regular gasoline. The charge rate and severity of catalytic-reformer operation 
was designated. The compositions of products including lead level in the motor gasolines were 
specified. 

In addition, shadow prices for each equation in the model and the relative unprofitability 
of each activity not in the basis was obtained. These data give valuable insight into the work- 
ings of the model and point out particular restraints which are costly. They also point up the 
relative merits of various activities which could be brought into the basis without probable 
serious consequences. 

Several variations in the solution of the problem have been obtained by parametric 
programming and by varying widely the values of the right hand quantities of the equations. 
Other variations in the solution have been obtained by blocking out selected activities. Addi- 
tional activities have been tested by the relative profitability criteria to determine if they had 
been originally built into the model whether the basic solution would have been changed. 

It should be understood that the coefficients associated with the activities, particularly 
the objective form, do not in any sense represent real and true values as used in the solution 
of the original problem. They have been deliberately altered in this presentation but are 
sufficiently similar to the original values that they should not distract from the principles 
involved in formulating the model. 


CONCLUSION 

To date Richfield has relied entirely on outside service centers to do the actual solution 
of linear programming models. The number and size of problems have not justified obtaining 
the large digital computers required for the solution to this type problem. It seems likely that 
Richfield and other smaller companies will be forced to continue on this basis for some time. 
This makes it possible for the people trained in formulation activities to concentrate their efforts 
on the intricacies of problem formulation and thereby hopefully develop more meaningful and 
valuable models. It is based on this background and the probable continuance of future efforts along 
these same lines which prompts the following statements on the current status of computing 
technology: 

1. There is a big need for speedy low-cost solutions. The linear programming models 
which are developing are steadily increasing in size. The end to this growth in size is not 
apparent. As the models get larger, the cost of calculating them increases very rapidly. 

There is a need for frequent solutions and, almost always, many re-solutions of the same 
basic problem. The re-solutions are required with only minor changes in coefficients or 
constants. These changes result from improvements in the validity of the model itself 
required to simulate actual conditions. There must be continual adjustments for changes in 
business environmental conditions as well as for the incorporation of steadily improving 
technology in the science of the business. 

2. There is need for complete printouts from the computers at the optimal solution. 
This, of necessity, requires not only the activities and magnitude of the activities in the basic 
solution, but also a shadow price for each of the equations (simplex multipliers) and relative 
profitabilities of the activities excluded from the basic solution. 

3. There is a need for flexible codes or computer programs so that better advantage 
can be taken of known feasible solutions, many of which, based on past experiences, may be 
very near to the optimal. Since variables with upper bounds occur frequently in problems, 
the inclusion of the upper bounding technique in the codes would be desirable to eliminate the 
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limiting equations. Curtaining techniques should be developed to facilitate exclusion of certain 
activities from the solution until an optimal is obtained with a printout of information at that 
point, followed by removal of the curtain to include all activities with the final optimal solution 
printed out. It should be possible to measure quantitatively the magnitude of the effect of 
removing one or more activities from a basic solution. While the relative profitability coef- 
ficient gives an indication of the magnitude of change which will occur by the removal of one 
activity and the substitution of another, it is often desirable to obtain additional insight into the 
model by determining the actual contribution of each specific activity. This is particularly 
true in areas near the optimal. Parametric programming techniques are useful for this 
purpose but the speed and cost of obtaining successive parametric programming solutions is 
excessive. 

4. There is need for further simplified methods for presentation of data to the com- 
puter, not only for the original problem, but for minor and major changes in the problem. 

5. It is desirable to be able to proceed with a minimum actual knowledge of the com- 
puter machine operation and a minimum knowledge of the coding or programming of the data 
calculations by the machine. In other words, for some purposes it is preferred to have a 
superficial understanding of the computer and how it works, leaving for limited sized technical 
forces primarily the job of model formulation. To be most fruitful, the efforts of a small 
group should be directed towards model formulation based on their particular business and 
its problems, rather than utilizing too much of their time in learning the intricacies of the 
computer itself. 

6. Since the petroleum industry is a process industry and many aspects of it are 
based on chemical change, it is entirely characteristic that many functions which must be 
built into linear programming models are of a nonlinear but hopefully convex nature. Many 
of the current characteristics of linear approximation of curvilinear functions or multiple 
recipes of alternate operations to cover the complex relationships introduce a very large 
number of at least variables and oftentimes equations to the model. It would be desirable if 
those working in the applied mathematics side of this business would devote more effort to 
means of handling these nonlinear functions in a more straightforward and simplified manner 
within the confines of the linear programming format. 





ON THE DISTRIBUTION OF MANHOURS TO MEET 
SCHEDULED REQUIREMENTS 


Saul I. Gass 
International Business Machines Corporation 


Washington, D. C. 


An engineering research and manufacturing company has its work organized by 
projects. The project number will be denoted by i, (i= 1,2,..., m), and it is known that a 
total of A; manhours must be applied to project i during a given sequence of months in order 
to successfully complete the project. For a given month j, (j = 1,2,..., n), the project plan- 
ning group knows how many total manhours, Bi, are available for use by the projects active in 
month j. At least qd; manhours and no more than eij manhours must be expended on project i 
in month j. If we define Vij as the total number of manhours assigned to project i in month j, 
the above verbal description can be expressed mathematically by the following constraints: 


n 
Oe nA 


m 
(2) 2 Vij $ Bi 
(3) o < dij £ Vig $ ey - 


It is assumed that for any two projects i and k, project i has a higher priority in the assign- 
ment of available manhours if i<k. It is also assumed that for any two months, p and q of 
a project, month p has a higher priority if p <q. In other words, the desire is to satisfy the 
requirements of project i before any assignment is made to project k (i < k) by allocating the 
available manhours to month p before any manhours are assigned in month q (p <q). The 
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computational problem is then one of finding a set of Vij that satisfy (1), (2) and (3) in a manner 
that also satisfies the priority restrictions. ! 

In devising a suitable computational scheme we are faced with two interrelated diffi- 
culties: The one being that for a given set of data (A,, Bi, di; ei) we have no assurance that a 
set of Vij can be found that satisfy the constraints; and the second being that even if the prob- 
lem is feasible, a set of Vij that satisfy the priority restrictions to any "suitable degree" may 
not exist. It should be mentioned that there is the implied restriction of not having, for a given 
i, the Vij fluctuating markedly between successive months. In what follows we will briefly 
discuss these aspects of the problem and describe a simple computational procedure that 
yields answers acceptable to the project planning staff. 

There are certain basic relations that must be satisfied if any feasible solution to con- 
ditions (1), (2) and (3) exists. Since ¥ij 2 qj; we must have 


and since 


j=1 i=1 


r mse 
cz &. 
a j-1 3 





The original mathematical statement of the problem was, find a set of yij that satisfy the 
constraints 


n 
a yij =Aj 


m 
Bj < Lyij s Bj 
i= 


re) es S tce Se. 
eer ij * Vij *eij 
and minimizes 


jae- Eval 
max |B; - Yi; 
j ; i=l a 


where B; and B; are, respectively, the minimum and maximum number of manhours available 
in month’ j (the difference is mainly due to available overtime hours) and B; = 1/2 (B; + B;). 
The Bj represent "desired" number of hours to be used in month j. The Bj in inequalities & 
are, in general, equal to the corresponding By. 
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As pointed out by Schell in his paper on multidimensional transportation problems [3], we have 
the following additional necessary conditions. Define 


m;; = min (A,, B;, ei) : 


We see that each mij represents the maximum value possible for the corresponding Vij? i.e., 


2 


Since 


we have the conditions of 


Since the lower bound conditions must be satisfied, we can rewrite the contraints of 
the problem by letting 


ye *% 


as follows: 


(1') 


(2') 


(3') 


Except for a statement of the objective function, the conditions (1'), (2'), (3') define a problem 
that is identical to the scheduling problem considered by Dantzig [1]. Dantzig demonstrates 





2 
For the data being used, we have found mjj= eij 
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how the problem may be transformed into a basic transportation problem and discusses a 
computational procedure which assumes an explicit objective function and a known first fea- 
sible solution. The problem with an explicit objective function is also equivalent to the 
capacitated transportation problem as formulated by Ford and Fulkerson [2] and can be 
solved by their algorithm. 

In the development of a suitable computational procedure for the problem defined by 
(1"), (2"), (3"), a number of things external to the mathematical model had to be considered, 
i.e., the number and size of the problems to be computed each month and the time allotted for 
the computations. Each month the project planning group is required to develop manhour 
allocations for approximately 360 such problems which involve, on the average, 12 projects 
that extend for 15 months. (Such computations represent a major repetitive effort of 20 
people.) If we wanted to use either of the schemes proposed by Dantzig or Ford and Fulkerson, 
we would have to develop at most 360 sets of cost coefficients each month that would yield 
acceptable answers, i.e., answers that conform with the priority ratings. It was felt that this 
would be a difficult and costly task as a lot of experimentation on subjective weighting would 
be required. These factors caused us to look for a computational procedure that would be 
simple and fast to apply. ; 

Our attempts were then directed to the development of a noniterative procedure based 
on a valid interpretation of the priority restrictions. We wanted to devise a simple scheme 
that would yield a feasible first solution that would also be acceptable as a final solution. 
Also, if for a few problems the procedure did not yield an acceptable feasible solution or a 
feasible solution at all, we would want the resulting information to be of some positive value 


to the planning group. It is felt that the procedure described below meets the above require- 
ments. It has been tested on a number of problems and has yielded acceptable answers in all 
cases. 


To describe the procedure we will refer to Figure 1 which represents a problem with 
four projects that extend for four months. Here the upper bounds, aj, are written in the upper 
left-hand corner of each box. The projects are ordered by priority, i.e., Project 1 has a 
higher priority than Project 2, etc. 








Project 









































Figure 1 
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We determine an initial allocation of the available manhours (this allocation will in all 
probability not be feasible) by starting the allocation with Row 1, Column 1 and continuing the 
allocations along the first row. Hence we find X44 = min (), hy, a44) = the maximum number 
of manhours that can be allocated to Project 1 in month 1. Ina similar manner x,, = min (; - 
X41» Mgr Ayq), Xyg = min (ry - Xyy - Xyp, hg, ayg) and xy 4 = min (ry - X11 -Xj9-%13,h4, 1 4)- 
This type of assignment is consistent with the priority ratings in that Project 1 is allocated 
the maximum possible hours. We then continue with Project 2 and determine: 


Xpq = min (To, hy - X44, 493), 
Xoo = min (To - Xp, he - Xz, Ag9); 
Xp3 = min (To - Xp) - Xpo, hg - X13, 49g) and 


Xpq4 = min (To - Xq1 - Xgq - Xp3, Ng - Xp 4, Agq4)- 


In general we have 


/ j-1 i-1 


Ea Hye 


The Xj are evaluated in priority order, i.e. » X. for j=1, » N; XQ; Fo? eee ee 
Xj io j= 1, ,n. If after the Xj are evaluated by (4) — all 


then this set of Xij represents an acceptable feasible eenuiten to the problem and the com- 
putation is complete. In general we would expect some r,> O and we are then forced 
with the task of adjusting the Xij in a manner that will tend to modify the low-priority 
allocations and also yield a feasible solution. Although we have no guarantee that the adjust- 
ment scheme to be described below will yield a feasible solution to the general problem 
described by (1'), (2'), (3'), we have found that for actual and experimental data the procedure 
has developed feasible solutions in all cases. It is felt that the adjustment scheme will yield 
acceptable feasible solutions in the greater majority of cases and the partial solutions to the 
remaining cases would represent excellent initial solutions to be further refined by hand. 

We will illustrate this scheme using the (4 x 4) example, Figure 2, assuming the Xj 


have been determined and defining hy = = hy - ps » ™ >0. We look at the rj in order (starting 


with i= 1) and find the first r, >0. Let us assume that it is r, 4? 0. We then search the h, 
in reverse order (starting with j =n) and find the first h 4 = 0. Let it be hy = = 0. This implies 


that the lowest priority Xj that can possibly be increased is X45, i.e., X44 and x44 must be 
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X42 | *43 
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ho = hg 


























Figure 2 


equal to their respective upper bounds. We wish to increase X49 by some positive amount 6 
and we see that 6 is restricted by the constraints . 


(a) X4o + 6< a4 


(b) r4- 620. 


Since we have added @ to Column 2 we must subtract a @ to maintain the balance that 


4 
be Xjo = hy . We first attempt to subtract @ from X39, the next higher priority allocation 
i=] 


in Column 2. To maintain the row balance in Row 3 we add a 6 to the adjoining low-priority 

' 
allocation in that row, i.e., X39. Finally, we have to subtract a 6 from hg to keep the 
accounting for Column 3 in order. We then have the additional constraints that 6 must satisfy 


(c) X39 - 6@>0 
(d) Xgq + 9 < Ags 


(e) : 


as pictured in Figure 3. 
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hy hg - 6 hy 























Figure 3 
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The maximum @ that satisfies the constraints (a), (b), (c), (d), (e) is 


' ' 
9, = min (ago - X49, T 4» X39 233 - X3q» hg) - 


If 95 >0 we have adjusted the solution as indicated in Figure 3 and, of course, if 5 =r, 4 this 
adjusted solution is a final solution. If r, . = r, .” 9, > 0, we again try to finda snitiee value 
of 6 by adding or subtracting 6 as above except that we now add @ to X34 instead of X39 and 
subtract 6 from hy instead of “ We continue finding 6 o® in this manner until we have 
shifted the adding of the @ to the last column. After this ~~ if the residual r 4 is still 
positive, we shift the above process to the second row (i.e., subtracting 6 from X59 and 
adding @ to X59) and repeat the evaluation of 9,'S; stopping whenever the residual r 47 0. If 
the process reaches x,,4 without the residual r 4 equal to zero, we then shift to the third row 
and first column, i.e., add 6 to X44) subtract 6 from X94 and add @ to X39, and repeat the 
above along the third row and up the first column until either the residual r 47 =0or X14 is 
reached again. In the first case a feasible solution has been reached and in the second case 
our process has been completed without reaching a feasible solution. 

Below is a sample computation:> 


Initial Tableau: 





Project || Jan Feb Mar 
3 2 



































+ 
0 


(0, =0 since 1+ 6 <1) 


3¥For this example > Aj= = Bj. 
i 5 
J 
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(@,=1 since 1 - 6 20) 


First Adjustment Tableau 
2) 
0 


1 
1 














y) 
a 


3 
1 























The next two steps yield 6 oo” 0 and they are not shown. 
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(6, =1 since 1 - @ > 0) 
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Final Solution 
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EXPERIENCES WITH THE BID EVALUATION PROBLEM 
(ABSTRACT) 


Hans Bremer, William Hall, and Maxine Paulsen 
National Bureau of Standards 


During the three years in which the Army Quartermaster Corp bid-evaluation problem 
has been solved, various techniques have been developed to handle the more common types of 
conditions that bidders might impose in the statement of their bid. In doing any contract- 
award problem the main difficulty is in establishing the cost matrix corresponding to pur- 
chasing and bidding conditions. In most cases the bid conditions require the solution of more 
than one transportation problem. In the code for SEAC, use is made of a programing device 
which permits consideration of all bid conditions simultaneously. The paper then describes 
various types of bid conditions and shows how these may be treated with the SEAC code. 








WAREHOUSING AND DISTRIBUTION OF A SEASONAL PRODUCT! 


Wm. S. Jewell 
Massachusetts Institute of Technology 


PROBLEM STATEMENT 

The basic problem of this study was to determine «DEMAND 
an economic means of warehousing and distributing a 
particular seasonal product. The scope of the problem PRODUCTION 
may be seen in Figure 1, which shows the total demand 
for the product (at the distribution level), and the planned 
production as functions of time. Because of the extremely 
high peak demand, pre-season production must be stored 
in warehouses, and then reshipped to supplement the peak | 
season production. as 

In addition to the problem of timing shipments to 
the warehouses, management was also concerned with Figure 1 - Total demand and 
the efficient shipping of the product into the various aaa —— saaieieaaiaaisi 
customer areas. Distribution to a given area could be 
made either by truck or by railroad from several different 
plants, subject to the receiving limitations of the various customers. The various shipment 
rates, warehousing costs, and production differential costs were available; with considerable 
effort, simple forecasts of customer demand by area by month for the entire year were also 
obtained. The production schedule for the year under consideration was available, and the 
amount of left-over stock in each warehouse was determined. 

The problem statement: Given all of the operation costs, the production schedule, and 
the estimated demand in each customer area for each month, how should the warehousing and 
distribution of this product be planned in order to meet the severe seasonal demand at least 
total cost for the entire year ? 





~ 





TIME 
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PROBLEM FORMULATION 

Because the distribution system under consideration had several warehouses and many 
customer areas normally assigned to each plant, the many possible combinations led to a 
problem of very large size, especially when considering the varying situation during each of 
the twelve months. Certain simplifications were computationally necessary, and also realis- 
tically desirable, since the first results of the problem solution were to be used only as a 
planning aid. 

Examination of the problem showed three types of shipment which occurred either by 
truck or railroad: 

1. Plant to customer area, 

2. Plant to warehouse, 

3. Warehouse to customer area. 
These shipments occur in various combinations during different parts of the year as illustrated 
in Figure 2. 

Because of the basic similarity between this problem and the elementary Hitchcock- 
Koopmans transportation problem, it was decided to use this formulation, rather than a more 
general and complex linear programming statement. 

The assumption made in order to use transportation problem techniques was that a 
fixed warehouse could be assigned to each plant-customer combination. This choice gave a 
unique route and associated cost for the combination of type (2) and type (3) shipments. 
Because of the large costs involved in warehousing and redistribution, the "least-cost ware- 
house" for each plant-customer combination was picked. This simplification effectively looses 
the identity of the individual warehouse, and only considers two types of shipment from a given 
plant to a given customer: the "direct shipment," and the "indirect," or "delayed shipment" 
(Figure 3). 


> 


\: 


PRE - SEASON SEASON PEAK 


Figure 3 - Simplified dis- 
Figure 2 - Types of shipment tribution scheme 


Another advantage of this formulation is that the problem of individual warehouse 
capacity is avoided. This omission is not a serious one since no solution to date has over- 
loaded a particular warehouse. Also, for the product under consideration, additional ware- 
housing is easily obtained. 

By arranging the plant supplies and customer demands in time sequence, one obtains 
an expanded transportation problem matrix which contains submatrices of the two types of 
shipment (Figure 4). Note that for m plants, n customer areas, and T time periods, a 
shipment matrix of size mnT? is obtained which has (m + n) T - 1 shipments as solution 
elements. 
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The physical arrangement of the sub- : DE a aad 
matrices automatically satisfies the produc - 
tiontotal from each plant during each period; nant 
the customer demand during a given month 
may be satisfied directly from production 
during that month, or from delayed (ware- 
house) shipment which originated during 
previous months. 

The costs to be used depend upon the 
position of the submatrix. In position (A), 
for example, the submatrix costs are the 
total of direct shipping costs and plant- 
production differential costs during the third 
month. Submatrix (B), however, has cost 
elements due to placing production from the 
second month in a "least-cost warehouse" for two months, and then reshipping it to satisfy 
customer demand during the fourth month. Normally all costs below the diagonal matrices 
are infinite, unless the possibility of back-ordering is allowed. In this case, the cost of con- 
sumer or distributor good-will must be evaluated. 

Having stated the warehousing problem in the standard transportation-type framework, 
all of the usual computational techniques are applicable. To determine an optimal production 
schedule, production limitations are used as supply, with a fictitious slack customer introduced. 
Stock on hand at the beginning of the season, or surplus production to be stored for another year 
may easily be introduced into the model. Shipment size limitations, simple probabilistic situa- 
tions, the cost of changing levels of production, etc. can also be added to the model with little 
difficulty. 


CUSTOMER AREA DEMA ND 


1 


2 
SUPPLY 


3 


PRODUCTION PERIOD 


4 


Figure 4 - Warehousing and distribution 
problem expanded shipment matrix 


COMPUTER FORMULATION 

Coincident with the formulation of the problem and the gathering of the necessary data, 
J. B. Dennis of the Electrical Engineering Department, Massachusetts Institute of Technology, 
- programmed the solution of the transportation problem on M.LT.'s Whirlwind I Computer, 2 
using a modified version of the stepping-stone method. 

The program contains 3500 instructions, 1000 of which are the actual stepping-stone 
method. The capacity of the machine solution is presently: 


(1) m < 127 | m = no. of plants 


(2) m+n < 401 = no. of customers 


n 
(3) No. finite costs + {ms 9) < 10,240 J = no. of row segments 
of infinite cost 





2The Whirlwind I Computer is a 16 bit/word digital computer with 2,048 registers of core 
Storage, and 24,000 registers of drum storage. The average speed is 35,000 operations per 
Second, although multiplication takes only about 12 psec. 
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Of course, in the expanded formulation, each original dimension is multiplied by the 
number of months during the season to obtain m and n for the computer input. Thus, this 
program is capable, for example, of handling a warehousing problem with four plants, thirty- 
six customer areas, for ten time periods. 

The unusual features of this computer solution are: 

1. Logical path searches are simplified by a unique method of storing elements 
in the basis. 

2. The criteria for introducing a new element into the basis is to pick the first 
element possible; + this results in more iterations per solution, but less passes through 
the cost matrix. 

3. Infinite element storage is simplified by means of "J" assignments which 
label large segments of infinite costs (non-allowable solution positions). 

4. It is also possible to obtain the alternate elements of the solution, or to use 
previously read-in costs with new supply and demand data. 

At present, there is no means for using a previously generated trial solution. This is not a 
handicap for small problems, but would be a desirable feature for large problems of a 
repetitive nature. , 

Computer experience with this program at M.LT. is shown in Table 1. Program read- 
in takes one minute; costs, supply and demand (5-place accuracy) are read in at about forty 
items a second. 


TABLE 1 
Computational Experience with Transportation Problem Solution 
on M.LT.'s Whirlwind I Computer 





Finite Data Input Plus 


No. of No. of Solution Ti 
Size Iterations /|Passes a a 


(%) (Min.) 








9 x 69 80% 150 1.5 
30 x 51 60% 314 3.3 


60 x 291 27% 1800-2400 























PROBLEM SOLUTION 
After the necessary data was gathered, it was formulated as above in a matrix of size 
60 x 291, which was approximately 27 percent finite costs. Computer solution time was a 
little over half an hour. 
Some idea of the solution results may be obtained from the solution cost figures below 
(Table 2). Projected Actual Cost, to which all other costs are normalized, is a projection of 
the previous year's actual distribution costs. The breakdown of costs includes: 
1. Transportation costs—railroad and/or truck. 
2. Warehousing costs—holding and throughput charges. 





Other criteria are being investigated; it appears possible to make substantial reductions 
in solution time in this manner. 
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TABLE 2 
Solution Costs for Warehousing and Distribution Problem 





Linear Programming Solution Costs 





Planned Production Schedule |Optimum Production Schedule 
Projected Actual Costs 





Experienced-Cost| Least-Cost Least-Cost 
Shipping Shipping Shipping 








Transportation | 100 84 73 71 





Warehousing 100 





Differential 100 








Total 100 86 79 67 


























3. Differential costs—production differential costs at different plants. 
4. Total cost. 

The first two solutions shown were run with management's planned production schedule. 
Two sets of cost figures were used: (1) Experienced-cost shipping uses the truck-rail ratio of 
shipments as presently experienced; (2) Least-cost shipping indicates that the most economical 
mode of shipment was assumed. The figures show that if it were always possible to ship via 
the most economical carrier, warehousing would cost more, but there would be a decrease in 
total cost due to the large saving in transportation costs. 

The third solution was run with the plant production capacities in place of the planned 
production schedule, thus letting the solution determine the most effective timing of produc- 
tion. The least-cost method of shipping was assumed, so that this solution determines the 
greatest possible saving available to management by using linear programming techniques. 

Naturally, it cannot be expected that all of the 14 percent to 33 percent saving indicated 
in Table 2 would be realized in practice. The vagaries of seasonal forecasting will always 
require a cushion of immediately-available warehouse stock; employment policy may not 
permit a production schedule which closely follows the seasonal demand. 

However, these preliminary results do point out to management the possibility of 
potential savings by (1) re-evaluating their present method of production scheduling, and 
(2) using linear programming techniques to determine warehousing and distribution policies. 
Even a few percent saving in distribution costs will pay many times over the cost of the 
research involved. At present, the management concerned with this seasonal product is 
undertaking a large-scale data-collection program, so that these and other techniques may 

« be used on a regular basis for overall planning and control of distribution. 
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PROBLEM EXTENSIONS 


The more general warehousing and distribution problem in which the capacities of 
intermediate warehousing are critical can also be solved using simple transportation problem 
techniques. Research on these and other multistage allocation problems is continuing at the 
Massachusetts Institute of Technology. 
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1. INTRODUCTION 

In the past few years various authors have proposed different methods for solving the 
transportation problem (e.g., see [1]). M. Gerstenhaber [2] has proposed still another 
approach which is novel to linear programming. It is based on the following interpretation 
of the problem. 

In a certain economy there are m producers who supply n consumers with a homoge- 
neous product. The je consumer has an ordering policy which states that if the unit price, 
tip quoted to him by the jth producer relative to the other producers exceeds a certain 
positive amount, r, he will not purchase anything from that producer. However, if the price 
does not exceed the threshold r, he will purchase a certain proportion, Piys of his total 
requirements, b,, from the jth producer. That is, he will purchase Xj = Pij b; units of the 
product from that producer. 

As a possible ordering policy we might assume that the jt consumer buys in propor - 
tion to the relative prices, ti - min, tig quoted by the producers. Thus, taking 


q,; = max (0, r + min, ty - ti) , 


a possible ordering policy would be 


(1) Pij = 45/2, es] R 


Now suppose that each producer has the objective to sell all his stock. Then, if the 
demand on the jth producer exceeds his stock level, aj, he will increase his prices. But, if 
his share in the market is small, he will most likely reduce his prices in order to get a larger 
share. 

We may look at this in another way. Suppose cy is the largest price the jth producer 
can demand of the jee consumer. Then his actual price may be written ty = Ci -Ais where 
4; 20. To increase his share of the market he increases A, ; to decrease it, he decreases A,. 

Now, Gerstenhaber has shown that if 5 a, = 2, b,, and Pi belongs to a certain class 
of functions oF defined in (1) belongs to this class), then for every positive value of the 
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threshold, r, there exists a non-negative i; for each producer which guarantees he will sell 
all his stock. In other words, for every r > 0 we may obtain a feasible solution to the system 


™_ 


25 y= b; 


Xj 20. 


Furthermore, he also shows that as r approaches zero, the feasible solutions approach 
the optimal solution to the problem of minimizing 


” ee 
subject to (2). 

Computationally, the method of solution consists in determining the producer prices for 
some value of r. This problem has been approached by means of a variation of steepest 
descent based on certain intuition about how producers and consumers might react in the 
situation depicted above. When a near feasible solution is obtained, the value of r is decreased 
and the process repeated. This goes on until r becomes sufficiently small. 

Although convergence of the actual process used so far has not as yet been proved, 
several nontrivial problems have been solved by means of it. It was observed that approxi- 
mately the same number of iterations were required for each problem, regardless of size. 

In the particular cases tried the number of iterations required ranged between 60 and 80. This 
is certainly an indication that the technique may be particularly well suited for large trans- 
portation problems. 

The apparently rapid rate of convergence can probably be attributed to two things. The 
first is the cut-off property of the functions Pij- These automatically price out of the market 
those producers whose prices are too large. For small values of r, this property selects 
practically all the Xj which appear with non-zero weight in the minimizing solution. The 
second property concerns certain monotone characteristics of the functions used (only the Pj 
defined in (1) have been used thus far). This will be discussed in greater detail in Section 7 
below. 

Thus, it is the purpose of this paper to extend Gerstenhaber's fundamental results to 
general linear programs with the hope that similar computational success will be obtained. 
Since the computing procedure is still in the experimental stage, we will limit our remarks 
to only a few computational aspects. The details of computation will be left for a future paper. 


2. THE PRIMAL AND DUAL PROBLEMS 
Given a set of m-dimensional column vectors = (b;), (0 <j <n), with Q 2 0, and 
scalars q; (1 <j <n). The fundamental problem of linear programming may be stated as: 


Problem 1: 
Find a set of scalars 2 such that 
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4 9-4, 


z; 20 


f= zy q; Zz = minimum . 
The dual of problem 1 may be stated as 


Problem 2: 
Find an m-dimensional row vector y' = (y;) such that 


rQas4 
g=y' a, maximum . 
We will make use of the following fundamental 


THEOREM (duality): If problem 1 (problem 2) has feasible solutions 
and f(g) is bounded from below (above) then problem 2 (problem 1) has 
feasible solutions and g(f) is bounded from above (below). 


Moreover, in this case 
min f= maxg. 


If either problem has feasible solutions, and its linear form is unbounded in the 
direction sought, then the other problem has no feasible solutions. 


3. NORMAL FORM OF THE PRIMAL PROBLEM 

In what follows we will require, among other things, that the primal problem be feasible 
and have a bounded minimum. These requirements lead us to a special normalized form of the 
problem. 

We proceed first to insure that the primal problem has a bounded minimum. We see 
that if Problem 1 is feasible then a sufficient condition for it to have a bounded minimum is 
that for some i, bi > 0, 0 <j <n, because then each variable is bounded both from above and 
below. That is, 0 < zs bi /dy; , 

However, if the primal problem does not satisfy this "positivity" condition, we may 
force the condition to hold by introducing into the system a constraint which bounds the sum 
of the variables, Z . That is, we select a large positive number N and let 


2; z+P=N, 


where p, a non-negative slack variable, is introduced with zero "cost." Although there may 
be some difficulty in selecting this bound in theory (especially if the problem really has an 
unbounded minimum), in practice technological conditions place natural bounds on the various 
levels of activity. Thus, it is generally not difficult to select N. 
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Assume now that b, >0, O<j<n. We may use this fact to make all the other ele- 
ments of the vectors positive by adding to each element of a certain positive multiple of 
its first coordinate. That is, we perform this same linear transformation on each of the 
vectors Q; . If we denote the transformed vectors by Q = (b,,), we may take 


' 
Now make the substitution 
2 = x aA bi /2j by) ‘ 

The primal problem then requires finding non-negative x;'s which minimize 

= 

he ie 

subject to 

z= Pj xj = Po - 


P; = Qa; bi 9 ae = qj /2; bij ’ Py = Q,/>4 Bio i 


Finally, we need to insure the feasibility of the linear program. This may be done by 
introducing m "artificial'' unit vectors into the system with high "cost," M. How large M 
should be is a moot point considering only a priori information. However, taking M > n max,c; 
seems reasonable in the practical case. 

In summary we obtain 


THEOREM 1: Let P; = (a; )), O<j<n, ande, = (55) > 1<k<m, be m- 
dimensional column vectors such that aij > 0 for all i, j, where 


wi aya, sjsn, 


. 1, isk 
ik 0, izk. 
Then problem 1 may be imbedded in the following linear program: For given scalars c 


7 
1 <j <n, find scalars x; and Yj such that 


(5) a, c;x,+M 2, y, = minimum 


subject to 
2 x P; + 2, Y; ey = Py , 


(6) 


ie y,20, 
where M >> max, 1 
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It follows as a consequence of (4) that 
(7) z x+ 2, yal. 


In other words, the solution of the transformed primal problem lies in the simplex defined by 
X; 20, y;2 0, and (7). 

We may, without changing the nature of the problem, include (7) in the constraint set. 
If we do this then the dual of the normalized primal problem may be stated: Find a scalar Wo 
and m-dimensional row vector w' such that 


(8) w+ w' P, = maximum 


subject to 

Wo + w' P; < c; 
(9) 

Wo + w' e; s M. 


By construction, the problems defined by (5), (6) and (7), and (8) and (9), both possess 
feasible solutions. 


4. A CLASS OF PRIMAL-DUAL FUNCTIONS 

It is easily seen that an arbitrary constant may be added to the linear form f without 
changing the nature of Problem 1. In particular, let this constant be -u' Q,> where u' is an 
arbitrary m-dimensional row vector. Substituting for Q,> f then takes the form 


(10) f= 2, (dj -u' Q) g,. 


If u' is a solution to Problem 2, then d,; - u' Q, > 0, for all j. Further, if u' is a maximizing 
solution, then from the duality theorem, f = 0 is attained for any minimizing solution to Prob- 
lem 1. 


By virtue of (10), (5) may be rewritten as 


(11) F (cj - w' Pj - Wo) Xj + 24 (M - w' e;- W) >» 


where W,» W' constitute a feasible solution to (9). In view of the linear dependence present in 
(6) and (7), a feasible solution may always be obtained by selecting w' arbitrarily and then 
defining W, as 

(12) Wo = ming, (c, - w' P,, M-w' e}) - 


Now assume that x and y, are functions of w and a "threshold" number r (> 0) with 
the following properties: 


a.) * (w, r), Yj (w, r)>0 


b.) 2 % (w, r)+ 2; y; (w,r)=1 
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c.) x; (w, r) = 0 if cy; -w' Py - Wy 2T 
y; (w, r)=0 if M-w' e;-w2r, 


where w,, is taken as defined in (12). 
That functions having properties a.), b.) and c.) exist may easily be shown. Let 


Uj=1 + Wy - ( - Ww P)), Vj=r+Ww,-(M-w'e,), 
uy = max (0, U)), v; = max (0, V;), 
A=, ujy+ 2 vy. 
Then the functions 
(13) x; (w, r) = u,/A y; (w, r)=v,/A 


clearly satisfy a.), b.) and c.). 
The functions defined above define a mapping which carries the set of dual feasible 
solutions into the space of primal solutions. We show in the next section that the image of 


this mapping contains a feasible solution to the normalized primal problem for every positive 
value of r. 


5. EXISTENCE OF FEASIBLE SOLUTIONS 
Define the functions 


(14) Fi (w, r)= 2) aj; x; (w, r) + y, (w, r), l<i<m. 


On the assumption that the primal-dual functions are continuous in w, we will show that for 
every positive r there exists a w > 0 with 


(15) =, Ww,=(m-1) (+ s)/, 


such that 


(16) F, (w, r) = a; l<icm, 


oO > 
where 
ony 
max; (M - c;) ‘ 


(Note that 4 #0 since a,, > 0 for all i, j). 

The functions F; Ww, r) define a mapping, G, of the m-dimensional w-space into 
another m-dimensional space whose coordinates will be denoted by gj - Thus, in order to 
prove the existence of feasible solutions, it is sufficient to show that every point of the 
simplex So in the g-space defined by g; 2 0 and 
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(17) {q*4%,°* 


is covered by the image of the simplex S,, in w-space defined by w, > 0 and (15). 
We first prove the following 


LEMMA 1. If for some i, w, = 0, then F, (w, r) = 0. 


PROOF: From (14) it is sufficient to show that y; (w, r)=0 and xj (w,r)=0, 1 <j<n, 
whenever w; = 0. From c.) this will happen when 


(18) cj -w' P; -W,21, forall j 

(19) M - w' e, - w, 29 

and 

(20) W, = min, (M - w' e,) = M - max, w, . 


6=(r + s)/A 


and consider the set of numbers Oy» 1<k <m, such that 


Clearly, the W, satisfy (15). Then 


max, Wy - W'P; = max, (a, + 8) - i aig (a, + 6)= 


(21) 
max, Me ~ Me Hy Ay + Oar O Arrt M-c, 


for all j. 
From (21) it follows that 


cj ~w'P; > M - max, W, 
Thus, (18) and (20) are satisfied. By construction 


max, W,29>Tr. 





42 J. E. KELLEY, JR. 


Thus, substituting (20) in (19) shows that (19) is also satisfied. This completes the proof of 
the lemma. 

From what was said above, we can prove the existence of feasible solutions by proving 
the following 


THEOREM 2. The image of Sy under G covers S.. 

PROOF: Lemma 1 says essentially that under G the vertices of é., are carried onto 
the vertices of Ss, , the one-faces into the one-faces, and so on for each higher dimension. 
This is easily seen as follows. The «th vertex of a, is defined by the vector w with Wy = 
(m - 1) 6, w,=0, i4k. Then Lemma 1 and (17) say that g, = 1 and g; = 0, i#k. A simi- 
lar analysis obtains for higher dimensional faces except, at this point, we can only say that 
the mapping is into. 

We proceed by induction to show that the mapping of vertices onto vertices implies our 
theorem. 

Assume that the image of every k-dimensional face of Sy covers a k-dimensional face 
of S, . Since each face is topologically a sphere, the restriction of G to each face has degree 
+1, depending upon the orientation. , 

Now, a subset of these faces form the boundary of some (k + 1)-dimensional faces of 
sy and 8, respectively. Let Sy and s, be two such faces such that the image of Sy is 
contained in 8, From what was said above, the boundary of S,, maps onto the boundary of 
S,, with degree +1. However, if the image of Sy did not cover all of s,, then it could be 
retracted to the boundary by projecting it from a point not covered. We could then deform 
the mapping to a point by retracting the image to its center. The mapping of the boundary of 
Sy onto the boundary of s_ is thus deformed to a point. It must then have degree zero, a 
contradiction. Thus, the image of Sy covers S,. This completes the proof. 


6. DEVIATION FROM OPTIMALITY 
Assuming that for a given value of r a solution to (16) can be obtained, how far from 
the minimum of (5) is it? This question is answered by the following 


THEOREM 3. Let W, Fr constitute a solution to (16). Then 


(22) Wot WP, <2, cy x (WT) + MZ, y, (W,T)<T+W,+ WP 


°° 


Thus, as r approaches zero the corresponding feasible solution approaches the optimal solution 
of (5) and (6). 


PROOF: From c.) it follows on substitution of W, and w for Ww, and w, respectively, 
that (11) is bounded below by zero and above by r. By rearranging terms (22) follows. 


COROLLARY. w,,, w approach the optimal solutioa of (8) and 
(9) as r goes to zero. 


The fact that the convex of (9) has a finite number of vertices allows the somewhat 
stronger 
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THEOREM 4. There exists r,, (> 0) such that for every r, 
O<r<¥8,, the corresponding solution of (16) minimizes (5). 
PROOF: Let w,, w be a maximizing solution to (8) and (9), and let 


S = {i, j |w, + w'P;=¢;, Wo + w'e; = M}. 


Then it is clear that there exists r, such that for r (<r), @,, B;, i, je S, 0< a, <7, 0<8;< 
r, there is another feasible solution Wo W to (9) with the property that 


ww = a 
Wo + WP =; r+a,, jeS 


<ey-r j¢s 


—_ _ , 
,_= > P € 
Ww, + We; M r+B; ieS 


Rewriting obtains 
jeS 


ifs 


ie S 


igs. 
This last result, Theorem 2 and the duality Theorem thus imply the conclusion of Theorem 4. 


7. SOME COMPUTATIONAL QUESTIONS 

There are many methods discussed in the literature for solving systems such as (16). 
Which of these, if any, is "best" in this case is still an open question. However, there are 
certain aspects of the problem of which any "'good"' computing scheme should take advantage. 
In this section we will attempt to point some of them out. 

In the first place we envision an infinite iterative process for computing w for some 
value of r. Assuming F; is a continuous function of r and w, we start with an initial approxi- 
mation, w) for some "small" initial To: and compute 


6, = Fw ,r)-a,, l<ism. 


Then, presumably on the basis of the result 6,0), a new approximation, wt), is calculated 

and the process continued. When the 6;'s become "small" (e.g., ar oF <e), t. is decreased 
in value to form Ty: Starting now with the last computed approximation w, the process is 
repeated. This goes on until either the current value of r is less than a preassigned tolerance 
or the difference in the value of (5) for two successive solutions attains a desired degree of 
smallness. At this point we accept the current feasible solution as optimal. 
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The success of the computation, of course, hinges on the method used for varying w 
from step to step. This, in large measure, depends on the primal-dual functions used. It is 
not necessary to select functions which require that (15) be satisfied exactly. For example, 
the Wj in (13) may be replaced by Ww; + k for any value of k without changing the functional 
values. Thus, in these functions we could take y w= 0, obviating the possibly large values 
of (m - 1) (r + s)/A. 

Any monotone properties of the 6, can certainly be used to advantage. For example, 
let 6; be a strictly increasing function of Wi and strictly decreasing function of Wj » S84. 
Our objective is to find a root of the system 5 =0, 1<i<m. For any W nota root, some 
of the 5's will be negative since, from (17), x, 5, =0, for any w. Let 5, (w) <0 and assume 
that we can always find y, such that 5, (W,,..., Wye + My, +++) Wm) =0. Let 8 ={i |5, < o}. 
Then compute My for each k in S. The substitution wi Wi» ifs, Wj = Wj + Mi, ieS, then 
gives a better approximation to the desired root. Clearly, this process is convergent, since 
at each step each positive 5 is strictly decreased. Thus, if it could be shown that eventually 
6; < 0 for all i, then from the fact that 2, 6; = 0, we must have that 6; = 0 for all i. Acon- 
vergent process for finding w might be constructed along these lines. 

Functions having the above monotone properties are analogous, as Dantzig has pointed 
out to the author, to Leontief systems. In such a situation fairly rapid convergence should be 
expected for even crude methods of approximation. Gerstenhaber has constructed a class of 
functions with this property. However, they have not been tried in any computation as yet. 

It should be noted that condition c.) of Section 4 has the effect of selecting a subset of 
the possible activities which the primal-dual functions consider good feasibility candidates. 

A considerable amount of computational efficiency may be gained by taking advantage of this 
fact. If at each step of the computation the change in Wi is bounded, then one can calculate 
the minimum number of future steps required before a particular activity, not in the current 
class of feasibility candidates, can enter the class. Thus, such activities can be avoided in 
computing the 5's for this minimum number of steps. 

Usually the artificial activities may be ignored. They were introduced only to make 
the proof of Theorem 2 possible. Ordinarily, if a problem is feasible and r is small, they 
will not belong to the class of feasibility candidates. It is not difficult to determine when they 
are actually needed. 

Finally, we point out that some computational advantage may be obtained whenever 
small changes are made in the original problem for some reason. By starting with the solution 
to the old problem only a few iterations should suffice to find the solution to the new problem. 

In the same way, if one knows the optimal solution to the dual problem, then for small 


r, only a few iterations should be required to compute the optimal solution to the primal 
problem. 
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A PRIMAL-DUAL ALGORITHM FOR THE 
CAPACITATED HITCHCOCK PROBLEM 


L. R. Ford, Jr. and D. R. Fulkerson 
The RAND Corporation 


1. INTRODUCTION 

An alternative method, called the primal-dual algorithm, was recently proposed [1] 
for solving the general linear programming problem. The algorithm starts with a feasible 
solution to the dual problem if available (otherwise a solution to a pseudo-dual problem is 
used); associated with this is a restricted primal problem. The restricted primal is then 
optimized using the simplex method, giving rise to a set of multipliers that can be used to 
obtain an improved dual solution. The process is then repeated. This algorithm may be 
viewed as a generalization of a combinatorial method discovered by Kuhn [3] for the optimal 
assignment problem and later extended [2] to the Hitchcock transportation problem. The 
primal-dual method becomes especially appealing in transportation-type problems, since 
special computing schemes—e.g., the flow algorithm of [2]—can be used to optimize the 
restricted primal, and no appeal need be made to the more general simplex process. 

The purpose of this note is to show how the Hitchcock problem with capacity con- 
straints on routes can be handled by using the flow algorithm for the restricted primal. Just 
as in the case of the uncapacitated Hitchcock problem, the primal-dual method for the capaci- 
tated problem appears to require not nearly as many iterations as does the straight simplex 
computation. This, coupled with the simplicity of the flow algorithm for the restricted primal, 
indicates that the computation proposed here is extremely efficient. 


2. THE PRIMAL AND DUAL PROBLEMS 
The primal problem is to 


. : : m n 
(1) minimize 2 a di; X54 


subject to the constraints 


(2a) 
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(2c) O< Xi < Cij P 


where a;, by» Ci and di are given nonnegative integers. We shall assume that 


which is clearly a necessary condition for the constraint set defined by (2) to be non-empty, 
i.e., for the problem to be feasible. This condition is of course not sufficient, because of the 
presence of upper bounds on the Xj . We know of no simple feasibility test for the problem. ! 
However, the algorithm will discover infeasibility if it exists. 

The dual problem to (1) and (2) is to 


m n m n 
(3) maximize a a; a; + a bj 8, + 4 a Ciy Vij 


subject to 


(4a) di, ~ @y - By - 7432 0 


(4b) <0. 


Yj 
Notice that a feasible solution to the dual is immediately available, e.g., a feasible solution is 
given by a, = 8) =7,,=0. 

In the sequel we shall sometimes denote the sets of indices i=1,..., m and j=1,...,n 
by M and N, respectively, and the set of pairs ij by MN. Subsets of M and N will be denoted 
by I and J, respectively, and the product notation used for these also. For example, IN will 
mean the set of pairs ij for ie I, j«€ N. Complements will be denoted by barring the appro- 
priate symbol; e.g., I will designate the complement of I in N. 


3. THE ALGORITHM 


We first give a description of the algorithm in the large, deferring discussion of the 
subroutine of Step 2 below to later sections. 
To start out, select the following feasible solution to the dual problem: 


a,= min dij ’ 


j 
?——o di; - @;), 


¥4y=0. 





The problem is feasible ifand only if the value ofa maximal flow inthe associated network 
(see Figure 1), where the ith source arc has capacity a;, the j sink archas capacity bj, and the 
intermediate arc ij has capacity Cij is equal to K. 
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Step 1. In terms of the dual solution for the cycle (the output of Step 3, below, of the 
preceding cycle), define the sets 


(6a) A={ije MN| y <ol, 
Let C=AUB. 
Step 2. Use the subroutine of the following section to solve the restricted primal 


problem: 


} 
(7) maximize a Xj 


subject to 


(8a) LX, say (ie M), 
(8b) GeN), 


(8c) $ X, < (ij¢ MN), 
(8d) (ijeC), 
(8e) %j = Ci (ije A). 


The output of this subroutine will consist of a matrix X = (x; ) which solves (7) and (8), 
together with subsets I and J respectively. If ps =. = K, nsssnsiallien X is a minimizing 


solution for (1) and (2). Otherwise, proceed to bs 3. 
Step 3. Define new dual variables by 


Ge), 
Ge), 


GJ), 
Ged), 


(ije IJ nC), 
(ije IN n A), 
(otherwise) , 


k = min Fave (di; -a@ i 783 ~7 4g) min in lrg ]- 
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Each of the bracketed quantities in (10) is to be interpreted as + if its minimizing set is 

empty. If k =, terminate; no feasible solution to (2) exists. Otherwise, repeat the cycle, 
' ’ ’ 

starting with Step 1 and using the a, By» Yi ‘ 


4. THE SUBROUTINE 

The subroutine for solving (7) and (8) may begin with any X satisfying (8). An efficient 
starting X is that generated on the preceding cycle. Corresponding to certain of the rows 
i=1,...,m (columns j=1,..., n) of X we shall define "labels" Aj (45) as follows. 

For those i such that )° Xj < a;, define 4, =0. Next, select a labeled i and scan for 


all (unlabeled) j's such that 
(11) ijeB and %j < Cay 3 


for these j's, define ,=i. Repeat until the previously labeled i's are exhausted. Then 
select a labeled j and scan for unlabeled i's such that 


(12) ijé B and %, > 0; 
for these i's, define A; =]. Repeat until the previously labeled j's are exhausted. Continue 


this labeling process, alternating between the rows and columns, until either 


(a) My has been defined for some j with )' Xj < b, , or 
j 


(b) no further labeling is possible. 
In case (a), alternately add and subtract h from the sequence 2 


13 _ g cee, ? 
(13) “uy * yh’ Ady “ile Ath’ “de » 


where AQ =0, 


=X = = = gia 
jy M;’ iy My Jo MY ip M iQ’ ? 


=. = 
jy, iy i, Mik ’ 


h= min [by - 2 xp “igh ~ ugh? “ugly? “yd ~ “yd,” 


Mil 7? hey’ FAs]? 


thereby obtaining a new solution to (8) that increases the form (7) by h. Then repeat the 
labeling process. Eventually case (b) will occur. We then let I be the index set of the 
labeled rows of X, J the index set of the labeled columns. 





2In words: Proceed, in the column just labeled, to the position indicated by its label and 
add h; then go, in this row, to the position indicated by its label and subtract h, etc., stopping 
only when a row labeled 0 has been reached. 





PRIMAL-DUAL ALGORITHM FOR CAPACITATED HITCHCOCK PROBLEM 


5. PROOFS 

The proof that the subroutine maximizes (7) for a feasible constraint set (8) may be 
found in [2], provided one makes the preliminary observation that (7) and (8) represent a 
maximal-flow problem in the directed network shown in Figure 1, where the capacity of the 
iM source arc is as ae , of the je? sink arcis b,=b.- Lc, , and the 

iT jsijea 4 j isijeA 

capacities of the intermediate arcs P,Q, are c;, if ijeB, zero otherwise. We also note one 
other fact from [2], namely, that the maximal-flow value in the network of Figure 1, i.e., 


for X; in the constraint set (8), is equal to the minimal cut value, which is given in terms of 
the labeling I and J by 


’ : b 
(15) z + ; 


Source 


In the following lemmas unprimed symbols refer toa given cycle q of the algorithm, 
primed symbols to cycle q+ 1. The matrix X= (x; ) and the index sets I, J are the outputs 
of the subroutine for cycle q; thus we are assuming that (8) represents a feasible constraint 
set for this cycle. Notice that this will certainly be the case for q= 1, as then A is empty. 


LEMMA 1. If ije BN TJ, then Xj = 0. 
PROOF: This follows from the fact that the relations Xj > 0 and je J would imply by 
(12) that ieI. 


LEMMA 2. If ije BN IJ, then Xj = Cys. 
PROOF: The proof is similar to that of Lemma 1, using (11) instead of (12). 
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LEMMA 3. If 0,8, 74; satisfy (4), so do a, By Yj 
PROOF: Table 1 accounts for all cases. 





i: 


TABLE 1 
Cases Occurring in the Proof of Lemma 3 





Index|,, _, | “iy ~ 71-9] - Yi 
ij ij 
Sets - d;; - @, - By - 74) 








Cals 0 
Cais 0 
Cals -k 
C ars k 





ANlJ 
Anis 
ANIS 
Ants 





BAlJ 
BNI 
BONIS 
BIJ 

















We see from the table that with k determined by (10), the dual relations (4) remain satisfied. 


LEMMA 4. If ije A', then Xj = Cay. 
PROOF: If ij¢A, this follows from the fact that the subroutine leaves X; unchanged. 
On the other hand, the table shows that any ije A' that is not in A is in BIJ, and the con- 


clusion follows from Lemma 2. 


LEMMA 5. If ijeC', then x,,=0. 
PROOF: If ijeC, the subroutine does not change Xij° If ij €C, it follows from the 
table that ijeB NIJ, and Lemma 1 applies. 


LEMMA 6. The constraint set (8) is again feasible for cycle q + 1. 
Moreover, X may be taken as a starting point for the subroutine of cycle 
q+1l. 
PROOF: It suffices to show that Xj = Cij for ije A’, Xj = 0 for ijeC'. These results 
are expressed in Lemmas 4 and 5, respectively. 








LEMMA 7. The dual form (3) increases by k(K - V- 1 C4;) >0 in 
passing from cycle q to q+ 1. A 
PROOF: It is clear from (6) and (10) that k> 0. If either k= oor L %j = K, the 
MN 








algorithm terminates. We are assuming, therefore, that k <a and )> Xij < K. 
MN 
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Now, from (9), 


a aj (a; 7 a) - b (6; - Bj) + Ci (Yj = ¥4) 


L 
MN 
pe 


. 2. ies a 2 
ImnA Y iste i] 
Using the max-flow min-cut theorem, we have 


V=La,-_Z ss + : ss 
1 * INNA mtia “ij * whie* 


Ly cy+ 2 
MINA IJnA 


shows that the left side of (17) is the same as the bracketed expression of (16). 
Since k(K - V- Dic.,) 21 if XY x,, < K, the algorithm terminates in a finite number 
A i MN 4 


of steps with either k= © or ) x,,=K. It now follows from the duality theorem that the 
MN 


1) 
original problem is infeasible in the first instance. That X is a minimizing solution in the 
latter case can also be seen from the duality theorem as follows. Since we have - Xij =a; 


and )° x,, 


~ ij = bj , and the relation Vij < 0 implies Xi = Ci , we have 


Fy" 8- se 
ae ele is Be 


But qj; -a; -B; Vi > 0 implies Xj = 0 (observe that iy < 0, qi - a, -B; iy > 0 is 


impossible in the algorithm), and hence the right side of (18) vanishes. Thus X is a mini- 
mizing solution. 
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TRANSLATING THE METHOD OF REDUCED MATRICES TO MACHINES! 


Bernard A. Galler and Paul S. Dwyer 
Department of Mathematics, 
University of Michigan 


1. INTRODUCTION 

The method of reduced matrices was designed to provide an efficient method of solution 
for general transportation (or distribution) problems of the Hitchcock type [1] in which unit 
parcels at origins are to be transported to certain destinations, subject to a specified matrix 
of costs for moving one parcel from each origin to each destination, so as to minimize the 
total cost. Several other problems are mathematically equivalent to this transportation prob- 
lem if one includes maximization problems as well as minimization problems. (Any maxi- 
mization problem can be reduced to a corresponding minimization problem by multiplication 
of the elements of the cost matrix by -1.) Maximization problems include the assignment 
problem (which now has the connotation of assigning specific individuals to specific jobs), and 
the personnel classification problem in which individuals or groups of individuals are assigned 
in groups to jobs. Other equivalent mathematical problems are available in the literature. 

Since the matrices of the grouped problems may be expanded to give square matrices 
with unit frequencies in each row and column, the essential mathematical problem may be 
interpreted as finding a permutation set of elements (one and only one element from each row 
and from each column) having optimal sum. In a grouped problem the grouping of this per- 
mutation set may be said to be a grouped permutation set. 

In the general transportation problem, as we use the term, we do not limit the discus- 
sion to two dimensions (as shown by the rows and columns of the specified matrix), but we 
consider a supermatrix (which we call a matrix) in k dimensions which has the form of a 
mathematical cell having equal numbers of rows, columns, layers, etc. The objective is to 
find the permutation set of these elements (resulting from the selection of one and only one 
element from each row, column, layer, etc.) having optimal sum. As in the common case with 
k= 2, grouping may frequently be applied to the permutation set as well as to the rows, 





l this report was supportedinwhole or inpart bythe United States Air Force under Contract 
No. AF41(657)-9 monitored by the Operator Laboratory Field Unit No. 1, Air Force Personnel 
and Training Research Center, Lackland Air Force Base, Texas. This is not an official publi- 
cation under the contract. Views or opinions expressed or implied hereinare not to be construed 
as necessarily reflecting the views or indorsement of the Department of the Air Force or of the 
Air Research and Development Command. Requests for additional copies should be addressed 
to: Engineering Research Institute, University of Michigan, Ann Arbor, Michigan. 
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columns, layers, etc., but alternatively the matrix of each grouped problem may be interpreted 
as a larger matrix utilizing unit frequencies. In this sense the mathematical problem for this 
general case, as well as for the special case with k = 2, may be said to be the determination of 
a permutation set having optimal sum. Since maximization problems may be reduced directly 
to minimization problems as mentioned earlier, the rest of this discussion is limited to mini- 
mization problems. 

This paper is intended to show how the method of reduced matrices has been translated 
into a workable computer program. Some of the details of the basic theory underlying this 
method have of necessity been omitted, and others have been modified in the process of adapt- 
ing the method to machine computation. In particular, the transformation described in Section 
6 (which is the most important feature of this method) is defined here in a way most suitable 
for automatic computation. 


2. THE OBJECTIVE OF THE METHOD OF REDUCED MATRICES 

The objective of the simplex method and other interchange methods as used in solving 
transportation-type problems with k = 2 is the determination of a feasible solution and the 
improvement of it by successive changes until the optimal solution is reached. The objective 
of this method of reduced matrices is to transform the matrix by operations which preserve 
the differences between the sums of the permutation sets. The resulting matrix consists of 
elements which are zero or positive, and in which some or all of the zero elements can be 
identified as those of a permutation set. Since the sum of the elements of this permutation set 
is zero, there can be no permutation set having a smaller sum. The objective of this method 
does not include the use of feasible solutions as does the simplex method, nor even the tentative 
partial assignment of frequencies to some of the elements presumably to be included in the 
solution as in Kuhn's method [11] for the assignment problem, or the Ford-Fulkerson* gen- 
eralization of Kuhn's method for the transportation problem, but concentrates on the deter- 
mination of transformations on the matrix which introduce additional zero terms. The aim 
here is to resolve the algebraic inconsistencies which appear in the application of the frequency 
conditions for each row, column, layer, etc. to the terms which constitute the trial permutation 
set. For many smaller problems to be worked by hand, especially with k = 2, many of these 
inconsistencies can be determined readily by inspection and the whole reduction is accom- 
plished informally. 

An objective of the method, then, is to provide techniques which are especially designed 
to handle this class of problems (determination of permutation sets having optimal sums) and 
thus avoid the artificiality introduced when these problems are adjusted so as to utilize meth- 
ods designed for more general problems. The artificiality of treating simpler cases as 
degenerate cases so that the formal requirements of the simplex method may hold, is a case 
in point. Also, the determination of a feasible solution (as has been pointed out by Schell [9]) 
becomes more difficult with increasing k. These difficulties are avoided by directing the 
solution toward the basic mathematical problem, the determination of a permutation set 
having minimum sum. 

The resulting solution is applicable to a wide class of problems in linear programming. 
An additional feature is that the machine process is the same for k = 2, 3, 4,... . No special 





*Ford, L. R., Jr., and Fulkerson, D. R., "A simple algorithm for finding maximal network 
flows and an applicationto the Hitchcock problem," The Rand Corp., Paper P-743, December 29, 
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program is needed for the two-dimensional case, although additional techniques are available 
for two-dimensional problems to be worked less formally by hand. 

In order to simplify the presentation, we are concentrating on a discussion of the case 
of k= 3. In general, the treatment of the k = 3 case exhibits the chief difficulties (the elimi- 
nation of negative solutions and the elimination of fractional solutions) of the general case. 

The k= 2 case is somewhat special in that there is no necessity for discussing the elimination 
of fractional solutions. 


3. THE GENERAL PROBLEM 

Let Cian be the known cost of transporting a unit item at origin i to destination h 
through intermediate point j. Denote the number of parcels at i by fi, the number to be 
delivered to h by fy and the capacity of the intermediate point by fi . It is understood that 
the capacities of the intermediate points are adequate, i.e., 


Ry Re Re 


Lft=L f= L t,=N 
int ? p21 J nos 2 


(3.1) 


where Ny, No, and Ng are the respective numbers of origins, intermediate points, and desti- 
nations. If we let Xian be the number of units, zero or positive, which are assigned to route 
ijh, we wish to minimize the transportation cost 


3.2 = . , 
(3.2) T itn Xijh Cijh 


(summation is over the entire range of the subscript unless otherwise indicated). These values 


of Xiih must also satisfy the marginal conditions (which we call the equations) 


a ijn = fj 


> =f 
ih *ijh ~ *j 


A Xijn = fh - 


If f; = fj =f, = 1 with Nn; = Ng = ng = N, the values of Xiih? as well as those of Cis form an 
Nx NxXN matrix. The set S of triples (i,j,h) with Xin = 1, form a permutation set. Then 
(3.2) becomes 
(3.4) y = by Cc: 

S ijh 
80 the problem is to find the permutation set S whose associated Cian have minimum sum. 


Any problem with non-unit frequencies can be interpreted as a permutation set prob- 
lem associated with an Nx Nx N matrix. The grouped results of this expanded problem may 
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be said to be a grouped permutation set consisting of the elements (i,j,h) with the associated 
non-zero frequencies Xijh . The problem is then the determination of integral Xijh > 0 which 
satisfy (3.3) and minimize (3.4). These values constitute a solution to the problem. The actual 
minimal sum may be obtained by substituting this grouped permutation set in (3.2). 


4. THE REDUCTION OF THE MATRIX 

Theoretically the type of transformation is simple. We subtract some constant from 
all elements of some particular row, column, or layer. Since in the ungrouped problem there 
must be one and only one selection from this row, column, or layer, the solution is not changed, 
Of course, the value of T is changed by the amount subtracted. This concept can be imme- 
diately extended to grouped matrices since the grouped matrix may be rewritten as an 
N XN*XN matrix with individual frequencies. Again the solution is not changed but the value 
of T is reduced by the total amount subtracted. Thus if a; is subtracted from each element 
of row i, which has frequency fj, the value of T is reduced by fa; . This process may be 
applied simultaneously to as many rows, columns, and layers as desired. 


5. THE MARGINAL TRANSFORMATIONS 

The first steps in this method of reduced matrices utilize transformations which result 
in non-negative matrices with at least one zero in every row, column, and layer. This can be 
accomplished by subtracting consecutively the smallest element in each row from all the 
elements of that row, the smallest element of each column from all of the elements of that 
column, and similarly with each layer. The resulting matrix is said to be reduced. 

In general there is nothing unique about a reduced matrix since the reductions in a k- 
dimensional problem can be performed in k! different orders, and the resulting reduced 
matrix may be different for each order. When using hand methods, it seems better to compute 
the amount of the reduction of T (the sum of the constants with frequencies taken into account) 
and to select the transformation having maximum reduction, i.e., first rows, then columns, 
then layers, or first columns, then layers, then rows, etc. For machine problems the sub- 
tractions in the order layers, columns, rows somewhat simplifies the process. After at most 
k steps, the matrix is reduced. 


6. THE GENERAL TRANSFORMATION 

Since this method of reduced matrices differs from other related methods most 
noticeably in the procedure used after the marginal transformations are completed, we are 
somewhat more formal in this section. Also, we state the problem in a slightly different way: 


Given a set of equations 


(6.1) Ax=b 


where A= a;,) is a p Xq matrix with p< q (here p= ny + No +... + my, and q= Ny Nog. . - Ny), 
find that solution x = (x, goons X,) in which each Xj is a non-negative integer and such that 


(0) 


q 
Ts c 
Ra Xj Cj 
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is minimized, where the Cj (0) are bounded. (In the transportation problem at hand, we have 
p+0, and a,,= 0 or 1 [cf (3.3)].) 

After the marginal transformations described in Section 5, the numbers Cc; (0) have 
peen changed into non-negative numbers c;, with a reduction in T of some seabdila To i. 


2 tst.08 aan 2 ae 
(6.2) oe % i . = i i’ 


and it is sufficient to minimize T - T ~ We assume that the c, are not all zero, otherwise 
any set of q positive integers Kyo ++ X constitutes a solution. 
Let Z= {j| c= 0}. We look for a solution of (6.1) of the form x= (x, ..., Xq) 20 

= 0 whenever j ¢ Z. (It is clear that for sucha solution, T - *, is minimized.) Let B be 
the matrix of coefficients of the system (6.1) subject to the condition x, = 0 for j ¢ Z. In other 
words, B contains the j-th column of the matrix A if and only if c,=0. By a succession of 
row transformations, B may be reduced to a matrix D = (d;;) (which is not necessarily unique) 
in echelon form, i.e., D has the form 


(* Dy 
6.3) De 


0 0 


where r < p, I. isthe rxr identity matrix, and D, may be empty. Then D=7B, where 
T= ;)) is the result of applying the same row —_— (in the same order) to the 


pX p identity matrix I: Let b' = 7b, i.e., “2, mb yr i=1,...,p. Then for a solution 


of the kind we wish, we have 
(6.4) Dx=7Bx=Tb=b. 


Suppose now that some row of D, say the i,-th row, contains only zeros, while bj # 0. 
o 


(This implies that the i,-th equation in the system (6.4) is not consistent with the other equa- 

tions under our restriction on the form of the solution. Such a row is called an inconsistency 

row, and bj is called the amount of the inconsistency.) By performing additional row trans- 
oO 





formations, if necessary, we may assume that for every inconsistency row i, bj > 0. We may 
also assume that among such rows, bj is maximal. 
fe) 


The existence of such inconsistency rows shows that a solution of the desired form is 
impossible unless the set 4 is enlarged. In this case, we define 6 by the equation 


" 
(6.5) 6 = min ; 
j¢Z ; ind 


r 
iy 
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p 
where % 4" 2 "ik aD j=1,...,q. Such a number 6 > 0 exists, since not all Cj are 


od 
zero, and if r; 43 0 for all j, we would have the contradiction (for any solution x = &,, eal 
oO’ 


%,) > 0 to (6.1)): 


q — __ P li 
02 2 ios) 2 hey 21 gk 9) > 24 “se *%, vi 


It can be shown that 6 is bounded from below by a number which depends only on the size of 
the original problem, since at each stage in the process the non-zero c; are bounded from 
below and the rj j are bounded from above by such numbers. 

oO’ 


We now define the transformation 
; (= C, - “ & Se 
(6.6) = ~ 6 “-i* j » 4 »4 
We observe that for j « Z we have cj “ee 0, since then r; a d; j7 0 by our choice of iy: 
oO’ to) 


Cc 


j : 
1 , and we have Cj = Cy -6r 


It is clear, also, from (6.5) that for some j, ¢ Z,6= 4° 
hh 1 41 ot 
Oo’ 


c; - c; = 0, so that Z is strictly increased to Z' by the transformation. Moreover, since 
1 1 


c, 20, j=1, 2,...,q, it follows that c; 2 0, j=1,2,...,q. If, for some j « Z, we had had 
j 


qd; <0 (cf Section 7 below) then c; =c, - 06 ry 4#- 6d, ; > 0, although c, = 0, so that in this 
oJ j od ij i 


case Z' would have fewer elements. 
Since bj can also be shown to be bounded away from zero by a number depending only 
fe) 


on the size of the original problem, it follows that T' = x. cj < T - T,, for any solution 
j= 
X= &, wise X,) of (6.1), since 


1 (°; : "A * 


= (T -T,) -@ 


ees Ag Te, * Fe: 


=1 


Let B' be the matrix associated with the system (6.1) relative to Z' as B was relative 
to Z. We may assume that the columns of B' consist of the columns of B with any new col- 
umns adjoined at the right. The matrix B' may be reduced by row transformations to an 
echelon matrix D' of the form (6.3), thus generating a transformation matrix 1’, so that 
D' = ' B'. We may also assume that in this reduction process those transfor mations which 
reduced B to D are applied first and in the same order. We shall show that D' has at least 
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one less inconsistency row than D, so that after a finite number of such transformations there 
will no longer be any inconsistency rows, i.e., the system of equations (6.1) will have a solution 
X= , ..., X,) such that X, = 0 whenever Cj #0, where Ci,---, Cy are the transforms of 
Cyr eres c_. It is clear from our assumption on the generation of 7' that all zero rows of D' 
are also zero rows of D, so that all inconsistency rows of D' are among those of D. We shall 
argue that the i,-th row of D' is not an inconsistency row in D'. 

Suppose D' = (aj) has . © = 0 for all j « Z'. Then the elements of the i,-th row of 


7B' can be expressed as linear combinations of the other rows of 7 B’, i.e., there are con- 


p p 
stants q,,...,q, suchthat )} 7,,a,= 5 YP aq, wy a, for all je Z'. Let X= 
. 4s m1 {of NS i7i mer tN 
&,, eee x.) be any solution of the system of equations obtained by dropping the inconsistent 
equations from the system B'x=b', so that, in particular, x; = 0 for j ¢ Z'. Then 


p 
(1'b); = b! - me q. b! = m TT. ~ 2 q; 7; 
40 to 44, tt mt igk Pk iZi, ey ik Pe 


r y 1 a, x 3 } > a 
= . X; - . 7 


p P , 
“ij orn Fs igk Aid” 24 re i Ti *3| 4 


= 0, so the i,-th row of D' is not an inconsistency row in D’. 


This process always transforms at least one inconsistency row to a consistency row at 
each step. (It may transform more than one.) We see then that a set of equations with m 
inconsistencies is reduced to a consistent set in at most m steps. 


7. ELIMINATION OF NEGATIVE SOLUTIONS 

When all inconsistency rows have been eliminated, a solution of the desired form can 
be found, and, in fact, is already at hand in the transform of the vector b by the row trans- 
formations used in the reduction described in Section 6. Although this furnishes a solution to 
the system (6.1), it is an algebraic solution, in that the x might be negative, and (except for 
k = 2) might be fractional. 

Suppose, first, that the i oth row of D contains at least one positive element and no 
negative elements, and + < 0 (we use the same notation as in Section 6, but we assume now 


that there are no inconsistency rows). By an additional row transformation, the i-th row may 


have all of its signs changed, so that now we would have bj > 0. We then perform the same 
fe) 


transformation on the elements cj as in Section 6, using the coefficients 7, ol as before. It is 


easily verified that the following results are obtained: 


427791 O-57-5 
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(1) For at least one j, ¢ Z, . =0. 


(2) Hd, 470 for some j « Z, then j ¢ Z' 


(3) All c; > 0 as before. 


(4) T - T, is decreased as before by an amount 6 b; which is bounded 
O 


from below. 
Condition (4) guarantees that the process will stop after a finite number of steps, 
resulting in a positive (not necessarily integral) solution, if whenever bj <0 (and the i,-th 
row is not an inconsistency row), all the elements of the i,-th row are non-negative. 


8. PERTURBATIONS 

One way of providing that the non-negativity condition of the previous section is satis- 
fied is to arrange to have D, be empty in (6.3). To insure this it is sufficient to insure that 
B have no more columns than rows, and this can be accomplished by insuring that all the ©) 
are distinct, so that there are exactly p-k zeros after the marginal transformations are com- 
pleted, and also so that each transformation introduces exactly one new zero. Another way of 
arranging to have D, empty in (6.3) is to temporarily replace each c; = 0 whose associated 
column falls in D, by a small positive quantity ej: Once a solution is obtained to the resulting 
system, the columns of D, may be returned to the system, and the corresponding Xj may play 
the role of parameters in the general solution. 

The latter treatment is effectively carried out in hand computations, while the former 
is quite easily accomplished by a computer at the start of the computation. In order to make 
the c (0) distinct, we perturb them by adding a (different) random number to each of them. 
These random numbers are added after first multiplying them by a sufficiently large negative 
power of two so that no accumulation of these numbers during the course of the computation 
can amount to as much as 1/2 of the unit in which the cj are expressed. The probability that 
the random numbers used have the same sum for two permutation sets having the same optimal 
sum is extremely remote. If necessary, however, another set of random numbers may be used. 

Such a procedure also modifies the problem so that it has a unique solution. By dropping 
these small random numbers we can then incorporate into the solution as parameters additional 
xj corresponding to the zeros thus obtained. 


9. THE ADJUSTMENT OF FRACTIONAL SOLUTIONS 

In a sense the attainment of a positive solution, fractional or integral, indicates the 
completion of the method of reduced matrices, even though an acceptable solution must be 
integral, since a matrix with a permutation set of zeros with fractional frequencies, subject 
only to the conditions stated earlier, can be reduced no further. This is shown by the fact that, 
if t is a common denominator of all the fractions of the permutation set of zeros, the corre- 
sponding cost matrix is a completely reduced matrix for a related problem with identical cj 
values but with frequencies tf. 

In a sense then we are faced with a second (combinatorial) problem. Using direct 
combinatorial methods we might select the elements corresponding to the zero terms of this 
matrix which is reduced as far as is possible, and then add in turn the terms corresponding to 
the smallest non-zero elements of this matrix and repeat the reduction process until a solution 
is reached. This technique is satisfactory for many practical problems since in many problems 
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one of these small elements does complete the solution. However it is unsatisfactory theoreti- 
cally, and it is easy to construct problems for which the proposal is obviously inefficient. 


The method developed here uses the matrix 7 resulting from the method of reduced 
matrices even though no further reductions of the matrix are possible. The elimination of 
fractional solutions by the addition to the solution of a single new element, as described below, 
has been incorporated into the program now in use on IBM 704. 


The method uses the particular fractional solution one gets by setting all parameters 
to zero in the solution obtained by the method of reduced matrices and modifies it into an 
integral solution with the addition of a new element which plays the role of a parameter. This 
means that the inclusion of the column associated with this element in the reduction of B to 
D leads to a reduced equation with zero in the column corresponding to the element, zero in 
the frequency column, and zero in the column corresponding to every element of the fractional 
solution. Furthermore, the coefficient of this element in the solution must be fractional and 
the fraction must be complementary to those of the particular fractional solution. These 
requirements give sets of necessary conditions which can be written in terms of the rows of 7 
and which must be satisfied if the element is to play the role of a parameter. Subsets of ele- 
ments of the matrix are determined from these conditions, the element (or elements) to be 
adjoined to the solution set is then determined from these subsets, and the general solution is 
obtained. 

Although the elimination of fractional solutions is an important feature of the method of 
reduced matrices, we have placed the emphasis in this paper on the determination of a positive 
solution, which may be fractional. Though improvements may be made in the method, we feel 
that the present techniques provide a practical and efficient method of solving this general 
transportation problem and that this method is generally superior, both by hand and by machine, 
to the simplex method even when k is as small as 2. 


10. AN ILLUSTRATION 

To illustrate the method of reduced matrices, we consider a k = 3, 4x 3x 5 problem 
proposed by Shell [9]. The matrix of cost elements Ciih is presented in Figure 1, with the 
associated row frequencies fi, column frequencies f,, and layer frequencies fh in the margins. 


Thus the element C394 = 34 has row frequency 27, column frequency 47, and layer frequency 
14. 





fj 16 





f 
i 
fh 14 17 17 27 45 14 17 45 14 





(2) 3 43 47 14 | 44 47 19 | 18 
6 17 49 15 20 | 47 1 9 | 22 
17 36 47 30 21 | 34 31 14 | 47 
32 16 12 17 15 | 47 14 27 | 16 




















Figure 1 - Reduction in T = 2(14) + 1(17) + 1(45) = 90 
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Although the process to be described below is almost the same as that used in the 
computer program written for the IBM 704, the perturbation of the cost elements described 
in Section 8 above is omitted (there is a unique solution to this problem and the perturbation 
process is not needed) so the process can be seen more clearly. 





16 47 





17 27 45 14 17 17 27 45 








42 47 13 | 42 47 (4) 23 18 
48 15 19 | 45 1 28 36 8 
46 30 20 | 32 31 17 2 13 
11 17 14 | 45 14 5 1 26 




















Figure 2 - Reduction in T = 1(47) = 47 


Figure 2 shows the result of subtracting the smallest element of each layer from each element 
of that layer (the elements transformed to zero are circled in each case). The numbers thus 
subtracted were 2, 0, 1, 0, 1, respectively, for a total reduction in T of 90 units (cf Section 4). 
A similar transformation performed on the three columns with frequency headings 16, 47 and 
57 yields a reduction in T of 47, since the only subtraction needed to obtain a zero in each 
column is a subtraction of 1 from each element of the second column. At this point the matrix 
is reduced (cf Section 5), and row subtractions are unnecessary. The result of the column 
subtractions is shown in Figure 3. 





16 





14 17 17 27 45 14 17 45 








o (3) 42 47 13 | 41 46 17 
4 17 48 15 19 | 44 0 7 
15 36 46 30 20 | 31 30 12 
30 16 11 17 14 | 44 13 0 25 




















Figure 3 - Reduction in T = 3(2) =6 


Now we pass to the more general type of transformation based on the elimination of 
inconsistent equations (cf Section 6). In Figure 4 the system Bx = b is presented, with zeros 
suppressed, the location of the c, = 0 in the matrix at the head of each column, and the equa- 
tions (3.4) identified at the left. (In the left-most column the first digit is the dimension and 
the second digit is the particular class within the dimension, so that the equation for Column 3 


is listed as 23.) The column headed by 112 is the result of the transformation about to be 
described. 
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111 123 222 232 234 333 424 435 112 






































Figure 4 - The system Bx= b 


It is not always necessary to reduce B to the matrix D in echelon form to find incon- 
sistencies. It is obvious that Eq. (21) and (31) cannot both be satisfied. (Since such obvious 
inconsistencies occur quite often, a search for the simple type that occurs here has been 
included in the 704 program.) In this case we can bypass the entire reduction process and go 
immediately into the transformation. The coefficients 1; oj 
(0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0). (This corresponds tothe eighth rowof the 12 x 12 identity 
matrix after the transformations 31* = 31 - 21 and 31** = -31* are applied to it.) By means 


needed to define @ in this case are 


c 
of (6.5) we find 6 = 112 _ 3 from Figure 3, and the amount of the reduction in T by this 
1 


transformation is 3(2) = 6, since bi = 2 inthis case. The result of this transformation 


fe) 
(which causes the column headed 112 to be adjoined to B in Figure 4) is shown in Figure 5. 





16 47 57 





17 27 45 14 17 17 27 45 17 





39 44 10 | 44 46 O 22 17 
45 12 16 | 47 O 27 35 7 
43 27 17 | 34 30 16 (4) 12 
8 14 11 | 47 13 4 O 25 




















Figure 5 - Reduction in T = 1(38) = 38 


Figure 6a shows the matrix D obtained at the next stage, except that the rows have not 
been rearranged to form an identity matrix as in (6.3). (Leaving the rows in their original 
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order facilitates the machine computation somewhat.) It can be seen that there are two incon- 
sistency rows, 33 and 35, with equal amounts of inconsistency. It does not matter which we 
choose to eliminate, and we simply choose the first, 33. The row of 1 (cf Figure 6b) corre- 


c 
sponding to this is also designated 33; from (6.5) we have 6 = 324 . 1, and in the resulting 
1 


transformation T is reduced by 1(38) = 38. The matrix one obtains at this stage is shown in 
Figure 7. 





111 123 222 232 234 333 424 435 b' 





14 
-12 
27 
31 





28 
27 




















Figure 6a - The matrix D 
































Figure 6b - The matrix 7 
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16 47 57 





17 27 45 14 17 17 27 45 





40 44 10 | 43 45 0 21 16 G) 45 
47 13 17 | 47 0 28 35 7 
44 27 17 | 33 29 16 0 11 0 
10 15 12 | 47 13 5 0 25 




















Figure 7 - Reduction in T = | (12) = 6 
2 


The process is repeated in Figure 8a, but this time there are no inconsistencies. (Note 
that two inconsistencies were eliminated with one transformation.) We have an algebraic solu- 
tion in the column headed b’, but there are four negative values. We use the largest (in size) 
to generate a transformation (cf Section 7). We apply one more row transformation, 12* = -12, 
to the matrix 7 shown in Figure 8b, and the second row, used to find @, is then 


(1, -1, 0, 0, -2, -1, 0, 1, 1, 0, 1, 0) ° 


c 
This time @ = _ = —, and the reduction in T is - (12) = 6. The result is shown in Figure 9. 





111 123 222 232 234 333 424 435 112 324 
































Figure 8a - The matrix D 
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Figure 8b - The matrix 7 





16 47 57 





14. 17 #17 #27) 45 14 17 17 #27 = «45 14 $17 #17 #27) 45 








0 0 40.5 44 10.5 | 42.5 445 0 20.5 16 17 0 445 40 20.5 
6 16 48.5 14 18.5 |47.5 .5 29 35.5 23 0 39.5 0 30.5 
15.5 33.5 45 27.5 18 33 29 #165 0 11.5 | 46.5 42.5 0 38.5 31 
31.5 14.5 11 15.5 13 47 13 5.5 0 25.5 | 16.5 16.5 24 2.5 0 




















Figure 9 - Reduction in T = 8(7) = 56 


Note that Co99 (which gave the solution value -12) is no longer zero after this transformation. 
It has been replaced by Ci39- We observe also that the number of negative solution values has 
been reduced by two. 


The process is repeated, as is shown in Figures 10a and 10b, with 6 = 
reduced by 8(7) = 56. The result is shown in Figure 11. 


111 123 232 234 333 424 435 112 324 132 
1 


= 8, and T 


C995 
1 



































Figure 10a - The matrix D 
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11 
12 
13 
14 





21 
22 
23 





31 
32 
33 
34 
35 























Figure 10b - The matrix 7 





16 


47 


57 





14 17 


17 «27 


45 14 17 #17 27 


45 


14 17 #17 


27 45 








0 0 
6 16 








15.5 33.5 45 
39.5 22.5 19 


40.5 44 
48.5 14 


27.5 10 
23.5 13 


2.5 
10.5 


42.5 44.5 0 20.5 
47.5 .5 29 35.5 
33 16.5 0 
55 13.5 8 


29 
21 





8 

0 

3.5 
25.5 





17 0 
23 0 
46.5 42.5 0 
24.5 24.5 32 


44.5 40 
39.5 0 


12.5 
22.5 
38.5 23 
10.5 0 





Figure 11 - Finally reduced matrix 


At this point the solution has been achieved, as is seen in Figure 12a, where b’ con- 
tains only non-negative integers, and all equations are consistent. 





111 123 232 234 333 435 112 324 132 225 





























Figure 12a - The matrix D 
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_ 


Nl] po] co 
Nie wl 


Ni— Nl 








1 
-1 -1 























Figure 12b - The matrix 7 


Figure 13 shows the solution values in the corresponding locations in the matrix (with 
zeros suppressed). The minimal value of T may now be computed in either of two ways. One 
may total up the successive reductions in T resulting from the transformations performed: 


90 + 47+ 6+ 38 + 6 + 56 = 243 units 


or one may apply (3.2) directly: 
14(2) + 15(2) + 2(0) + 2(0) + 2(1) + 38(1) + 2(3) + 25(2) + 13(2) + 7(9) = 243 units 





16 


47 


57 





14 17 17 27 45 


14 17 17 27 45 


14 17 17 27 45 





14 








2 





15 





13 
2 





Figure 13 - The solution 


The computation time for this problem on the IBM 704 was 1.3 minutes. 
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CONSTRAINT MATRICES OF TRANSPORTATION-TYPE PROBLEMS ! 
(ABSTRACT) 


I. Heller 
THE GEORGE WASHINGTON UNIVERSITY 
Logistics Research Project 


There are probably many ways of defining what is meant by transportation-type prob- 
lems. In common usage is the following definition: A problem is of transportation type if its 
constants and variables can formally be interpreted as the constants and variables of a trans- 
portation problem. 

The present note starts from the other end: A transportation problem is a particular 
linear programming problem; particular because of its special features regarding existence 
of integral-valued solutions and simplified computational methods; these features are reflected 
in particular properties of the constraint matrix which can be taken to define the term "trans- 
portation type."' A study of classes of matrices with such properties is the content of this 
note. 

The constraint matrix of the Hitchcock-Koopmans transportation model has the prop- 
erty that whenever a column is a linear combination of linearly independent columns, the 
coefficients are integers (and therefore necessarily +1 or 0). This property permits a sim- 
plification in the computational procedure of the "simplex method" as was shown by G. Dantzig 
[1], and assures that all basic solutions are integral valued whenever one integral-valued 
solution exists. C. Tompkins and the author [2] have shown the mentioned property to hold 
for an extended class of matrices. J. Kruskal and A. Hoffman [3] proved the property (which 
they termed "unimodular") to hold for certain classes of incidence matrices related to partly 
ordered sets. This speech before the Conference presents recent contributions to the subject 
comprising the following statements: 

» 2 by, bo, eared b, are linearly independent vectors in a linear vector space, then 

the set 


D={+b;, b; - by} (i 4 j; i, j=1, 2, — 


is unimodular. Or equivalently: 





1 Work done under the sponsorship of the Office of Naval Research. Reproduction in whole 
or in part is permitted for any purpose of the United States Government. 
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If the affine span of 


~~ * 


is of dimension n (that is if the end points of the b; are not in an affine subspace 
of dimension n-1), then the set 


D = {b, - by} Gi 44; i,j=0, 1,...,; n) 


is unimodular. Or equivalently: 

. The set D of edges of a simplex (each edge taken in both orientations and inter - 
preted as vector) is unimodular. 
D is maximal for its dimension, in the sense that the span of D does not contain 
a unimodular set properly containing D. 
Theorem. A unimodular set of dimension n contains at most n(n + 1) elements 
(not counting the nullvector). 
Theorem. A unimodular set of nfm + 1) vectors (not counting the nullvector) is 
of type D, that is, the set of edges of a simplex. 
For n > 4 there are n-dimensional maximal unimodular sets that are not of type 
D (they necessarily have less than n(n + 1) elements). 

For further details the reader is referred to the author's article [4] which contains 
also the proofs of above statements, while the following remarks intend to illustrate their 
significance. 

Since every subset of a unimodular set is unimodular, the interest lies on the deter- 
mination of all maximal unimodular sets. Note that a maximal unimodular set contains with 
each vector also its negative. 

Further, a maximal unimodular set A determines a certain class of maximal uni- 
modular sets, namely the class of all images of A under nonsingular linear transformations. 
Therefore, to describe all maximal unimodular sets amounts to finding one representative in 
each class. Clearly, for a given dimension n, the number of distinct classes is finite. 

Concerning these classes of maximal unimodular sets: Statements (1) and (2) confirm 
that the sets of type D constitute one class, which by (4) is the only class whose sets contain 
n(n + 1) vectors. Statement (5) confirms the existence of other classes, which by (3) consist 
of sets with less than n(n + 1) vectors. So far no attempt has been made to determine these 
latter classes. 

Geometrically, let A be a set of vectors in Euclidean space such that the span of A 
has dimension n and A does not contain the negative of any one of its vectors (note that a 
maximal unimodular set contains with each vector also its negative). Then one may consider 
the question whether the vectors of A can be so arranged as to constitute part of the edges of 
a simplex. An obvious necessary condition is that A be unimodular. Statement (4) says that 
this condition is also sufficient, provided A contains n(n + 1)/2 vectors; Statement (5) shows 
that it is not generally sufficient when A contains less vectors. 

For interpretation in a free Abelian Group see [5]. 

While all sets of class D can be obtained from any one representative by nonsingular 
linear transformations, those containing the unit vectors (and hence in coordinates having all 
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entries equal 0, + 1) admit a rather simple construction. A set U containing the unit vectors 
is obtained from an arbitrary set A of class D by choosing a basis in A and representing each 
vector of Ainthat basis. To choosea basis amounts to 
choosing a tree in the 1-dimensional graph of the 
simplex edges. All sets U are obtained by choosing 
all possible trees (with oriented and numbered 
segments). 

Figure lis an example ofa tree and the asso- 
ciated set U; the j-th column of U corresponds to the 
segment j; the negatives have been onitted. 

Figure 2 is a star whose segments are so 
oriented as to yield the set D of Statement (1) if the 
b; are unit vectors. 

Figure 3 is a chain, so oriented and numbered 
as to yield a set C which consists of all possible 
columns witha string of consecutive 1's (the negatives 
omitted). 

For computational aspects of findinga star or 
chain basis in a given unimodular set compare [6]. 











O 
) 
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Figure 3 
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ON THE ASSIGNMENT AND TRANSPORTATION PROBLEMS 
(ABSTRACT) 


James Munkres 
University of Michigan 


In this paper we presented an algorithm for the assignment problem which is a variant 
of H. W. Kuhn's so-called Hungarian method [1]. We also gave a generalization of it to the 
transportation problem. Since a detailed exposition will appear elsewhere [2], we shall 
content ourselves here with a few general remarks. 

The basic data of the problem is incorporated in a matrix called the rating matrix or 
the cost matrix. In the Hungarian method and its variants, the general procedure is to trans- 
form the matrix by adding and subtracting constants from rows and columns of the matrix, 
until a matrix is obtained which contains zero elements in a certain crucial arrangement. In 
the process of determining what numbers to add and subtract, one needs to distinguish certain 
rows and columns, which are usually called covered rows or columns. One also needs to 
distinguish certain zero elements, usually done by means of asterisks (and primes, as well, 
in the present algorithm). 

The solution procedure, involving many additions and subtractions, is often rather 
tedious if carried out by hand. The Operations Research group, ERI, University of Michigan, 
has constructed a special-purpose computer to aid in this procedure. The author's algorithms 
were developed for use with this computer. 

The computer consists of a 20-by-20 square array of teleregisters on a display panel, 
each one capable of displaying any integer from 0 to 10. There is a telephone dial attached 
which is used to set up the matrix on the panel initially. Next to each row and column are 
push-buttons—pressing one button will cause each number in the corresponding row (or column) 
to be diminished by 1; pressing the other will cause each number in the row to be increased 
by 1. 





Next to each teleregister are two lights—the white one is used to stand for an asterisk 
and the green one for a prime. They are controlled by individual switches. Finally, interlaced 
across the matrix vertically and horizontally are metal rods, painted black on one side and 
white on the other. One distinguishes a row or column by flipping the corresponding rod so 
that its white side is exposed. One could more easily distinguish a row or column by merely 
placing a distinguishing mark next to it, but the use of colored rods is better for our purpose, 
for in the course of solution it is desirable to be able to distinguish easily between elements 
which lie in distinguished rows or columns and those which do not. 
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The main purpose of the computer will be to serve as a demonstration device, to 
accompany lectures on the assignment and transportation problems and to illustrate various 
methods for their solution. It is too limited in size for many problems, since the largest 
problem it will handle is a 20 by 20 assignment problem, or a 9 by 19 transportation problem. 
One hopes, however, that the computer may provide valuable information relative to the feasi- 
bility of constructing a similar semi-automatic computer of a more useful size. 
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A QUADRATIC PROGRAMMING PROCEDURE! 


Clifford Hildreth 
Michigan State University 


1. THE INITIAL PROBLEM 
Consider a programming problem of the following form -- 


Problem I: max 6 (x) = x'Cx+ d'x 
x 


subject to Gx2h. 


x is an nxl (column) vector whose elements are real variables. C is a negative defi- 
nite matrix. G is an mxn matrix, d’ is 1 xn, h is mxl. The elements of C, G, d', h are 
known constants. If the underlying practical problem requires that some or all of the elements 
of x be non-negative, this can be incorporated in the above formulation by giving G an appro- 
priate identity submatrix and h a corresponding null subvector. 

Several algorithms have been suggested for problems such as I (see Frank and Wolfe 
[1] for an example and a short bibliography). All are iterative and some have the advantage 
of terminating in an exact solution in a finite number of steps. Compared with these, the 
procedures discussed below have the disadvantage that, except for special cases, the approach 
to a solution is asymptotic. However, this is at least partially offset by the ease of proceeding 
from one step to another, and Procedure I has seemed to work reasonably well in the several 
instances in which it has been tried. Variations in the procedure were suggested by a step in 
the proof of convergence for Procedure I and, to the writer's knowledge, they have not been 
applied. 


2, SOME EQUIVALENT PROBLEMS 
By a theorem of Kuhn and Tucker [4], a vector x is a solution to Problem I if and 
only if x and some z (of m components) solve the following minimax problem. 


Problem I: min max y (x, z) = x'Cx+ d'x - z' (Gx - h) 
z x 


subject to z20. 
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For any z, the maximizing x can be found by differentiation. 


(2.1) 


oF Laee d'-2z'G. 
ox 


Equating the partial derivatives to zero gives -- 


(2.2) xmax _ zc (G'z - d). 


Substituting this expression for x in (x, z) yields -- 
x 1 ,,' -1,' anol at oe | ' 
(2.3) vq, 8) = -2 GC G'z -2d'C*G'z+ d'C*d)+h'z. 


Let -- 


(2.4) = -7Gc7"G' ve sacha’ +h’, 


then, Problems I and II may be solved by solving Problem III below and applying (2.2). 


Problem III: min 9(z) = z'Az+ b'z 
Zz 


subject to z20. 


Two cases are distinguished. If G has full row rank, then A is positive definite. This 
case will be called Problem Ila. If G is not of full row rank then A is positive semidefinite 
with diagonal elements aj; > 0. This will be called Ib. 

We assert -- 


(Al) Problem Ia has a unique solution. 


Call the solution zmin Its existence is assured by the facts that 9(z) has a lower bound and 
the set of admissible z is closed. If there were distinct minimizing points, 
the convexity of 9(z) would imply -- 


(2.5) o(u) <o(zminy for any uso gin, y - ) gmtin with O<a<1l1, 
thus denying that gmin gmin are solutions. 


3. THE FIRST PROCEDURE 
Take an arbitrary z0) 20. Obtaina z (1) by minimizing 9(z) with respect to each 
component of z in turn while holding other components fixed at their last obtained values. 


z) is similarly obtained from z( , etc. In general, z(P+ 1) is found from z®) by the 
relations -- 
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z,°+ 1). max (0, w, P+ 1)) 


wi” + 1) is the value of z, that is obtained by setting pas = 0. In taking the partial 
Z:. 
i 
derivative and obtaining z+ 1) the preceding components are held fixed at their just- 


obtained values in z+ 1) The subsequent components are held fixed at their values in z®). 
For these fixed values of the other elements, 9(z) is minimized with respect to non-negative 


Z; by setting z,°* 1) = w,P* 1) if the latter is non-negative and z,°* 1) = @ if w,?* 1). 0. 


(3.1) and (3.2) define a continuous vector function which will be written -- 

(3.3) z+ 1) _ pz) y , 

Procedure I consists of choosing an admissible 20) and repeatedly applying the operator P. 
4. PROOF OF CONVERGENCE FOR PROBLEM Illa 


Successive applications of the operator P generate a sequence of vectors, {2®)}, and 
a sequence of scalars {o(z ). We wish to demonstrate -- 


(A2) If A in Problem I is positive definite, 2®)} _, gmin 


It is first useful to show -- 

(A3) P(z) = z=z= gmin | 
Take any 220, z7 zmin, then 9 (z™min -o(z)< 0. Let zmin _ 7 4 h and note -- 
(4.1) o(2™iM _ 2'Az+ 22'Ah+ h'Ah+ b'z+b'h, 


(4.2) 9 (2) = (2) + (22'A+ b')h+h'Ah, 
(4.3) 9 (2M) _ 9(z) = (2) h+h'Ah. 
Since the second term on the right of (4.3) cannot be negative, the other term must be 


negative, This implies that 9(z) can be reduced by changing one coordinate of z and the 
operator P will do this. 
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{2} is contained in the bounded set [z / 9(z) S o(2)) ] and therefore has a limit 


point. Let z” be a limit point and { alt]} be a subsequence approaching z”. To suppose 


ZZ z™iN can be shown to lead to a contradiction thus establishing (A2). From (A3), if 


af zmin then -- 


(4.4) o(z”) - o(P(z™))=e>0. 


The uniform continuity of 9(z) in the bounded region containing {z®)} assures that, for any 
such e there exists a 5 > 0 for which -- 


(4.5) ha-2* 1 < 5=>/ oz) - o(z*) /< e. 
From the continuity of P, there exists a p > 0 such that -- 

(4.6) Mau-2* [<< p=>/ P(z) - P(z*) /< 6 
and for any p, there exists a positive integer, R, such that -- 


(4.7) r>R—/f zit] _ M<p. 


7 
Take such an r, say r*, and let p* be the index of zit ] in the original sequence, {2)}. 
For such a p*, 


(4.8) / o(2®"*)) - (P(e) /<e 
and, therefore, 

(4.9) o(zP**1)) < 9¢z%) 

which contradicts the definition of {2} 


5. PROOF OF CONVERGENCE FOR PROBLEM IIb 

It will be recalled from Section 2 that Ib is the special case of Problem III in which A 
is semidefinite with aj; > 0. A problem of this form does not necessarily have a solution, but 
in many cases the investigator will know from the practical context that a solution does exist. 
If the investigator is uncertain, a procedure developed in the paper by Frank and Wolfe [1] 
may be used to test the existence of a solution. 

In general, solutions to problems with semidefinite quadratic parts are not unique. Let 


zMin be any solution and gmin _ » (zmin) be the least feasible value of the criterion, 9(z). The 


following assertion will be proved: 


(A4) If the matrix A of Problem III is positive semidefinite with positive diagonal 
elements and if the problem has a solution, then the sequence {9 (z®) )} generated 
by Procedure I converges to the minimum, 9 
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Note that the i of (A3) in Section 4 also applies to the semidefinite case. In these 
circumstances j9(z“’)} is monotone decreasing and bounded below (by gmin) hence it con- 
verges. Let » be the limit point. 


To suppose o” + gmin leads to a contradiction. Under this supposition we have -- 


(6.1) o® - gmin = >0. 


The function P(z) which is defined by (3.1) - (3.3) above may be regarded as a composition of 
functions of the form -- 


(6.2) P(z)= PL, Piya --+- Pp Py @) 


where P, P; (z) = P; (P; (z)) and P; (z) is the vector obtained from z by holding all components 


except z; fixed and substituting for z, the non-negative number which makes the value of the 
function » as small as possible. 
Let A, (z) be the amount the criterion is reduced by applying the operator P; to z,i.e., 


(5.3) A; (z) = 0(z) - o(P; (z)) 20. 


min 


Consider any solution gmin and any non-negative z. Set z = z+ h as in Section 4 and 


rewrite (4.3) -- 


(5.4) o(z) - 9 ™m . - (2) -h'Ah. 


Since h' Ah is non-negative 


(5.5) -(22)n = x ~ h; 2 (2) . g min 
i= 


Oz dz. 


1 


, dg 
Changing z; to - would reduce 9(z) by an amount equal to - 3z,°5,- By definition 4; (z) 
i 


is the maximum reduction in 9(z) attainable by changing z;, therefore 


a9 
(5.6) A; (z) 2 - 22, h; 


and 
(6.1) 4:2) 2 ofa) - o™™. 


The relation (5.7) is interesting partly because it can be used at any stage of the iter- 
ation to obtain a lower limit for 9™", specifically -- 





C. HILDRETH 
m P 
(48) ofa) - Ya, (a) s min 
l= 


for any element of {2}. It also follows from (5.7) that -- 


(A6) For any z 20, 3 k3 A, (z) zi (9 (z) - gmin) : 


In particular, under the supposition of (5.1), for any element of {z®)| there is an 
index k such that 4, (z®)) >= . If P, were applied to z) this would obviously deny (5.1). 
However, P. is applied after Pi, Po, wag Pyy and further examination of the matter is 


needed. For convenience, let 
*_ (p) 
(5.8) z=Py_y Pyig-:-: Po Py (z*’). 


Now for p sufficiently large, it will be shown below that 2”) and z* can be made arbitrarily 


close. A, (z) is continuous and for 4 z) _ 2* / sufficiently small, A, (2®)) and A, (z*) 


r 
differ by an arbitrarily small amount, say < — When this holds 9 (P+ 1)) <9 (z®) ) - _ 


But since this must hold for all p beyond a certain point, the convergence of jo”) )} is 
contradicted. 


It remains to show that / z) - z* / can be made arbitrarily small. If {o ))| is 
to converge, all A; (z) that actually enter must approach zero. In going from z®) to z+ 1) 
the operator P; is applied to a vector which agrees with z+ 1) in the first i-1 components 
and with z®) in the last m-i+ 1 elements. Call this vector % and let z= P, (2). 


A, @) = o(2) - o(P, @) ) 
(5.9) 


, 2 22 . as 
A, (2) = aj; @" - 2°) + 2(@; - 2) (2, ™5* +) 


From (3.2), using the notation of this section, 


1 1/2 oo 
(5.10) w,P* ye (2, ayy 2; + >. 


Combining (5.9) and (5.10) -- 
(6.11) A, (2) = aj; (& - %) (B; + %, - 2w, PP) ) , 


(5.12) A, (@) 2 a, (@, - 2)". 
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The inequality is justified by the fact that either Z; = w, P+ 1) 2 0 which makes the 
equality in (5.12) hold or Z; = 0, w, P+ 1). 0 which makes A, (z) greater than the magnitude 
on the right. / Z; ~ Z; / can be made less than any « > 0 by making p so large that 

2 
A (z) < Vai; Ee . 


6. POSSIBLE VARIANTS OF THE PROCEDURE 

The procedure described in Section 3 has been applied to several rather simple prac- 
tical problems. The author [3] has reported an example of Problem Ila and Freund [2] has 
solved an application of Problem IIIb to programming under uncertainty. These and a few 
unpublished examples do not give much basis for comparison with other procedures. It is not 
hard to imagine circumstances in which the number of iterations could be very large (about 
250 is the worst case the author has encountered) but, as noted before, each step is easily 
constructed. Of course, the investigator may sometimes be able to help by making a careful 
selection of 20). 

Assertion (A6) suggests a possible modification of the procedure. One might construct 
z P+ 1) . P,. (z®)) where k is chosen so that A, 2) 2 a, (@)) for i= 1, 2,...,m. That 


is, instead of applying the elementary one-component a P; in turn, one could look for 
), at each stage and apply it. This 
‘ min 
9(z) -9 
m 


the most effective one, that one which maximizes A i (z 


would insure that the criterion would be reduced by at least at each step. How- 


ever, looking for the most effective elementary operator is a somewhat larger task than per- 
forming an iteration with the original procedure. 

One could, of course, apply Procedure I at the beginning and use the modification along 
the way if convergence seemed unduly slow. Another modification that might be worth con- 
sidering would be to draw random numbers from 1 to m at each stage to determine the 
sequence in which the elementary operators would be applied. 
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COMPUTING PROCEDURES FOR PORTFOLIO SELECTION 
(ABSTRACT) 


Harry M. Markowitz 
RAND Corporation 


Let X; be the amount invested in the jth security. Let Ms be the expected return on 


the i » security; 0 iv? the variance of the jth security; and Oia the covariance between the jth 
and jt security. The expected return on the portfolio as a whole is 


.th 


The variance of return on the portfolio as a whole is 


n 
= “1 ms "ij X; x; . 


Portfolios are selected subject to constraints 
pm ayy X; = b, 
X% 20. 


(In the simplest case we have only one equation 1X, = 1.) The portfolio selection problem is 
to find those portfolios whose (E, V) combinations are efficient (see figure). 

The "critical-line" method for solving this problem was presented in the author's 
paper "The Optimization of a Quadratic Function Subject to Linear Constraints," Naval 
Research Logistics Quarterly, Vol. 3, No. 1 and 2, pp. 111-133. A more elementary expo- 
sition of the method will be presented in a forthcoming Cowles Foundation monograph. 

The "critical-line" computing procedure is similar to the simplex method for solving 
linear programming problems. A slight modification can convert a simplex machine code into 
a code for the critical-line method. 
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LOSS OF ACCURACY IN SIMPLEX COMPUTATIONS 


Walter Jacobs 
United States Air Force 


1. INTRODUCTION 

This paper deals with the solution of linear programming systems on digital computers 
using the simplex method. The number of decimal places or significant figures that can be 
carried along in such a computation is limited, and it is generally not possible to guarantee any 
specific degree of accuracy in the solution. 

In order to discuss how the loss of accuracy might be controlled, it is necessary to 
consider in detail how errors propagate. For this purpose it is useful to think of a hypothetical 
computer, performing the computation in parallel with the actual computer, but using exact 
values throughout. However, this parallelism may turn out to be impossible. If the compu- 
tation procedure has branch points, and if the choice of branch depends on the value of a 
quantity subject to round-off error, the approximate and the exact computation may follow 
different paths. Where the path chosen in the approximate computation is invalid, we say that 
a wrong decision has been made. Branch points where wrong decisions might occur we call 
critical steps of the computation. 

As long as no wrong decision has been made, the approximate and the exact computation 
differ only in the values that arise, and these differences are due solely to round-off errors, 
which behave according to well-known rules. When the approximate computation takes an 
invalid path, it is no longer possible to say in general what is happening to the accuracy of the 
values that are obtained. To see this, we will examine the inverse form of the simplex method 
in detail. No change in principle would be involved in similarly treating some other form of 
the simplex method. 


2. THE SIMPLEX METHOD — INVERSE FORM 
The algebraic statement of the form of linear programming system we are considering 
is the following: Minimize the expression 


(1) ao1 Xj + Ano Xo + --- + AnD X, 


subject to 
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aoy Xy + agg Xo + ...+ ay, X, = by 


Ag, Xy + Ago Xp + --- + Ag, X, = by 


and 

(3) x20, I) -— FF 

It is assumed that all b, 2 0, which can always be arranged by changing the signs of all coef- 
ficients in appropriate rows of (2). 


The computation procedure replaces the above statement by the equivalent problem: 
Maximize Xp subject to 


Xo + (@p, X1 + Ang Xy +--+ ap, X) =0 
Xio1 + @qy %1 + AyQ Xp +--+ + Ayn A) = PY 


X42 + @oy %1 + Agg Xp +--+ + Ann %) = Dg 


Tne * @u1 X1 + ane Xt: 


Hj 20 


+k 2 0, 


(5) 


where we define 


= - (by + by +...+b,) 
(6) 


yj = - Go + agp t+ --. + any), G2 3,8... ,@. 


The method consists of a succession of cycles, each of which begins with a tableau in 
the following form and ends with the tableau that starts the next cycle: 
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Index Value Inverse of Basis 





Yo 
‘4 


Vo 


where the variables x are those which occur in the basis, and all other variables x; are 
i 


zero. The entries in the index and value columns at the start of the first cycle are 
(7) = i , " (i= 1, 2,..., m) 


and Vo = 0. Furthermore, the starting form of the inverse matrix is the identity matrix of 
order (m + 1). 

The following are the steps in going from the tableau that starts a cycle to the tableau 
that starts the next cycle. For the moment, exact calculation is assumed, and no provision is 
made for degeneracy. 

I. Test Vy: If Vv, < 0, proceed to Ila. If v¥,* 0, proceed to Ib. (Step I will then be 
omitted from all succeeding cycles, which will start with Step Ib.) 
Ila. Compute q;, Gj=1, 2,...,n) as 


(8) dj = Bio 40; + Bij a4j + oer # Bim amnj ° 


Find s, the first index satisfying qd. = min d; . 
Test d.. If d. = 0, end the problem as infeasible. If d. < 0, proceed to Ila. 





. Compute Yj» (i= 0, 1, 2,...,m) where 


¥i = Bio 205 + Biz 415 + --- + Bim 4ms- 


Test each ¥j- If  : 0, ignore that value of i. If y; > 0, compute w; = v;/¥;- Find r, 
the first index satisfying Ww, = min w, (among those i for which w was computed). 
Proceed to IVa. 

. Compute the tableau for the new cycle as follows: 
The index column is unaffected, except for the entry i, in row r (r is the index of 
Step Ila), which is replaced by s (the index of Step Ila.) The other entries in that 
same row v,, B.,; are replaced by 





ri 
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' 
vp = V/A, 
By = By (k= 0, 1, 2,..., m). 
For all other rows i/r, (i=0,1, 2,..., m), 
Bi = By - 2 Bry (k= 0,1, 2,..., m) 
2 = ¥i/Vy - 


Having completed this, return to Step I with the new tableau. 
IIb. When the result of Step I is v,2 0, compute 


(12) ds = Boo 29; + --- + Bom mj » Gah 2... ,8. 


Find s, the first index satisfying d. = min 
Test d.. If d. = 0, end the problem as solved by the basic solution contained in the 
current tableau. If d. < 0, proceed to Ib. 

- Compute y, by (9) as in Step Ila. 
Test each y,. If y, <0, ignore that value of i. If y, > 0, compute w, = v,/¥;- 
Find r, the first index satisfying W, = min Wj: 
If there is no i with y; 0, end the problem as having no upper bound to Xo - In the 
contrary case, proceed to IVb. 

. Compute the tableau for the new cycle using (10), (11), and return to Ib with this 
tableau. 








3. THE CRITICAL STEPS 

The procedure just outlined contains two types of operations: First, computations such 
as (8) or (10), and second, decisions. These decisions with which we are concerned are the 
steps which have been underlined in the procedure described in the preceding section. 

If these underlined steps are reviewed, they will be found in every case to have two 
characteristics. First, they all involve a decisior about which of two or more branches will 
be followed in continuing the procedure. Second, this decision depends on comparisons 
involving computed quantities which are subject to round-off error when the procedure is 
being carried out using digital arithmetic.! 

The underlined steps are the critical steps of the computation, because at these points 
the calculation might result in choosing an invalid branch of the procedure. For example, 
round-off error in d, might result in continuing the procedure after the exact calculation has 
stopped at the optimum. In such a case, the basis ultimately obtained as the solution will 
almost certainly not be a correct one, and the difference between this and the exact solution 
can no longer be called "round-off error." 





lit is this last point which really distinguishes the designated steps from other ''transfers" 
such as occur on an internally sequenced computer when the last term of an expression like (8) 
has been computed 





ch 


(8) 
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We consider in detail the kinds of wrong decision that can arise at each of the under- 
lined steps in the procedure of Section 2. For this, we assume that no wrong decision has 
been made prior to the step under discussion. Furthermore, we suppose that there is avail- 
able an upper bound ¢ for the absolute value of the round-off error at the point.* We will use 
an asterisk to designate the approximate value corresponding to a quantity occurring in the 
exact solution. 

The first kind of wrong decision arises from an error in v, during Step I. If ¥%,* 
but vi < 0, we fail to recognize that a basic feasible solution has been obtained. In practice, 
most computer codes test for vi > -€. This slightly increases the risk that a feasible solu- 
tion will be assumed to have been reached when it actually hasn't, but a wrong decision of the 
latter type is usually much less serious than the other: 

Similar remarks apply to the test of d, = 0 in Step Ia or Tb. Again, the termination 
of the problem will be made to depend on d® > -€. Since this test and the previous one become 
crucial in every problem, all practical simplex codes make some such provision for round-off 
error. 

Turning to the remaining underlined steps, we note that one of them, the determination 
of s in Step Ila or Ib, does not properly rank as critical, because the validity of the compu- 
tation is not affected when an s* is chosen for which d.+* is not minimal, so long as d.* < 0. 
Next, the test of y; <0 in Mi or Ib will lead to an invalid path only if a wrong choice of r 
results; otherwise, the effect on the computation is momentary. Consequently, the choice of 
an incorrect r is the only type of wrong decision which remains to be considered. It is easy 
to see that the basis produced by such a wrong decision is actually infeasible. 

There are three different ways in which a wrong choice of basis can arise. With r* 
denoting the wrongly selected row of Steps Ila or IIb, we may have: 


(a) y's 
(b) y* 
r 


(c) ws 


The first of these is the most serious, since it completely destroys the validity of the compu- 
tation. It can often be avoided by trying a new value of s whenever _. < €, provided that 
r 


there is another s' with d,1< 0. Case (b) can only occur if there is some i for which 
0 > y; >-€ 


and at the same time 


* 
v, <€(1+ w",). 
i r 


This shows that (b), as well as (c), involves "near degeneracy." 





2We will not consider in this paper how to obtain such a bound 


427791 O-57-7 





94 W. JACOBS 


4. REDUCING THE LIKELIHOOD OF A WRONG BASIS 

If the bound € to the size of the round-off error is sufficiently small, the likelihood 
becomes negligible that a computed number less than € in absolute value is actually different 
from zero. Since a wrong choice of basis can be made only if this unlikely event occurs, we 
can eliminate wrong choices in the great majority of problems by suitably controlling round- 
off error. 

One way to do this is to use double-precision arithmetic throughout the computation. 
Judging by RAND's experience, this is successful although it may be costly in machine time. 
Another approach being investigated at the Air Force is to use single-precision arithmetic 
until one of the three cases of the preceding section comes up, and at that point reduce €, by 
improving the inverse (B;,), until the uncertainty about r disappears. The usual method of 
iteration to improve the inverse will also give an estimate of€ . We do not have any informa- 
tion yet on how this would compare in average machine time with the double-precision 
approach; in favorable cases it would seem to be faster. 

Instead of controlling wrong decisions by such methods, the problem might be cir- 
cumvented by using a solution procedure like Relaxation, which does not have invalid branches. 
Unfortunately, relaxation methods are still not competitive with the simplex procedure. It has 
been suggested that relaxation might be used to improve the inaccurate result obtained from a 
simplex computation, but I see no good reason to believe that a relaxation starting from a 
wrong basis would converge quickly. It appears that in linear programming systems of fair 
size, the only practical way at present to obtain answers of acceptable accuracy is to use a 
simplex method coded to avoid wrong decisions at the critical steps. 





SOME METHODS OF COMPUTATIONAL ATTACK ON 
PROGRAMMING PROBLEMS OTHER THAN THE SIMPLEX METHOD! 
(ABSTRACT) 


C. Tompkins 
University of California, Los Angeles 


This paper contains an exposition of methods originally due to Kaczmarz [1] and 
Brown [2] for the solution of problems in programming, including some problems in games. 

The methods are recommended for use only when the complexity of the system to be 
solved is great in terms of the versatility of the computing equipment available. This relative 
complexity may be due either to a particularly complex problem (with some nonlinear aspects, 
for example) or to the use of computing equipment of limited versatility. 

The Kaczmarz scheme for the solution of linear algebraic equations is one of projection 
orthogonally onto hyperplanes represented by the equations to be solved. The projection is 
usually made onto the hyperplanes in cyclic order, the point which is the result of the pro- 
jection at one stage being the point to be projected in the next. For extensions of this scheme 
to inequalities and to systems of higher algebraic order see [3], and for related work see the 
bibliography to [3]. 

The Brown scheme for the solution of two-person zero-sum games depends on each 
player in turn playing a best strategy against the average strategy used by his opponent in a 
history of moves. Robinson [4] proved that this method of fictitious play converges to good 
average strategies for both players. The convergence has been observed to be slow; some 
acceleration seems to be possible by redefining the game and limiting the fictitious plays to 
mixed strategies encompassing the strategies which seem likely to be the correct limiting 
strategies. 

The paper notes that the Kaczmarz scheme converges slowly for systems which are 
almost singular, and it also notes that the relations describing a game with the value filled in 
approximately is almost singular. The slow convergence is exploitable here to permit fast 
convergence to a value which approximates the solution of the game: that is, the Kaczmarz 
method for the noncompatible system of inequalities which result from game relations into 
which an approximate but incorrect value has been inserted, does lead to a series of 





The preparation of this paper was sponsored bythe Office of Naval Researchand the Office 
of Ordnance Research, U. S. Army. Reproduction in whole or in part is permitted for any pur- 
pose of the United States Government. 
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approximations to good strategies from which better approximations to the value may be com- 
puted. This scheme gives good convergence with modest computing demands. 

A particular example is cited in connection with problems of the type met in dynamic 
programming; this example is found in [5], [6] and [7]. It is in many ways similar to eigen- 
value problems, and the same technique could be applied to calculation of eigenvalues. 

The example consists of a game with a parameter in the matrix: the payoff. matrix has 
the form A= B=tC, where all the elements of the matrix C are negative. The game is to be 
solved under the additional condition that the value of the game is to be made zero. 

A modified Brown attack is recommended to solve this problem. Each player makes 
the most pessimistic estimate of the value of the parameter t which is consistent with his 
observations during the history of play, and under this assignment of value to t plays in 
accordance with the Brown method. In practice a few hundred plays seem to give enough 
information to permit use of the Kaczmarz method for acceleration and the more precise 
determination of strategies and the value of the parameter t. 

The paper will be submitted for publication at a later date. 
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EDITING LARGE LINEAR PROGRAMMING MATRICES 


P. M. Thompson 
Hanford Atomic Products Operation 
General Electric Company 


Large linear programming matrices may involve thousands of coefficients. It is a 
most difficult task to determine, record, and transcribe all of these coefficients into a machine 
medium without making an error. Yet the entire purpose of the matrix and all of the work 
behind it are served only to the extent of the accuracy of this work. The increase in reliability 
of results fully justifies all of the purposeful editing and proving that can be applied to the 
matrix. A secondary motivation is found in computing costs that are saved by eliminating 
computing reruns. 

There are three fully distinct steps in the process of arriving at a matrix of coefficients 
in machine language: (1) Representation of the real world by a mathematical (linear program- 


ming) model; (2) determination of the values of all of the coefficients of the matrix, and (3) 
transcription of the matrix coefficients into the machine medium. The model of a given real 
world situation will appropriately reflect its activities, commodity flows, goals and require- 
ments, physical and policy restraints, and the effective interrelationships. There must be 
value judgments on all of the activities. These are written into the objective or functional 
equation. 


OBJECTIVE EQUATION 

It is with respect to the value functional that some of the most important introspective 
questions must be asked, for it is this equation that will determine the form of the optimal 
solution of the matrix. The model can be no wiser than its functional. There are grave 
dangers that symptoms of a problem may be under analysis rather than the actual root prob- 
lem. Problem areas are first detected via their symptoms. It is imperative to delineate the 
entire problem area and properly locate the specific problem to be analyzed within the context 
of the entire area. This frequently involves going up the control hierarchy of the organization 
to levels higher than the level where the symptoms were observed. It is only through such a 
process that the optimization of the linear program model can be made completely compatible 
with the central objectives of the entire parent organization. 


ADEQUACY OF THE MODEL 
Deter mining the validity or adequacy of a mathematical model is a difficult process. 
At such time as one feels that the developmental work on a model is completed (i.e., the model 
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is finished) several things can be done. If time permits, the model can be put away for a 
period of time, then re-examined after some of the background thinking has faded. The effect 
of each equation should be determined as fully as possible working from the structure of the 
equation rather than following again the same chain of logic as used in originally writing the 
equation. The detailed operation of the model can be described orally to an uninitiated but 
critical mind. Defects will likely become apparent in this discussion. Study of optimal solu- 
tions is most important. Since the solution vector has many elements it is essential to array 
the solution in tableaus that reflect interrelationships of variables. A thorough detailed 
examination of early solutions is required, checking all variables for reasonableness and 
consistency. Defects in the model usually show up quickly in solutions. The more penetrating 
this analysis, the fewer solution runs will be required to evolve a satisfactory model. It may 
be worth while to examine the alternate vectors, P., (if any) that are eligible for the optimal 
solution (i.e., 6 °° P; = 0) to see if their inclusion would be reasonable from the physical point 
of view. 


GENERATION BY MACHINE 

Once the model is written, it should be studied to see how many of the elements can be 
generated by machine methods from the logical interrelationships between the rows and col- 
umns. Punched-card equipment can be quite effective in such generation. Both the location of 
the element within the matrix and the value of the element can often be generated. There will 
be cases where for some sections of the matrix, the i-j addresses of the elements can be 
generated even though the coefficients themselves cannot. In these instances, the i-j indices 


can be generated and listed for convenience in inserting the values of the coefficients manually 
for subsequent key punching. In some problems almost the entire matrix might be generated 
by machine. By proper application of standard machine techniques, the possibility of error 
can be virtually eliminated from this process. Machine generation should be given the greatest 
consideration because of its inherent accuracy. 


WRITING OUT THE MATRIX 

For those sections of the matrix that cannot be generated, it is important that they be 
written out on large sheets of cross-section paper, with the numerical values of the coefficients 
written carefully. The matrix can be broken into sections for convenience. A linear program- 
ming machine code must provide space for at least a five-digit name for each vector. Mnemonic 
names must be used suggesting the significance of each vector. These names along with row 
numbers and a notation of the physical significance of each row must be included on the matrix 
charts. In some cases it may be worth while writing out the sections that were generated 
mechanically for the array of the full matrix gives powerful insight into the logical structure 
of the matrix and the interrelationships of the variables. It is another check point of value on 
the construction of the model itself. 

Getting the numerical values of the coefficients on paper is the weakest link in the 
entire chain and from this point of view it is worth while writing out the entire matrix. Once 
the numerical matrix is written, it is possible to almost guarantee the transcription to machine 
language without error. 

After the coefficients have been calculated from the physical constants of the real world 
they should be thoroughly verified. They can be independently calculated by another person. If 
possible, they should be set aside to cool before recalculation. Where possible, a backwards 
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calculation should be made going from the coefficients to the original constants. The coeffi- 
cients should be checked by common sense. Do partial capacities look right? Do they add up 

to appropriate full capacities? Do rates and limits look reasonable? All possible cross checks 
should be made. 


CHECK SUMS 

It is advantageous to transcribe the matrix to punched cards, and then convert by 
machine from cards to the appropriate high-speed machine input medium. Only one matrix 
element should be put on each card. This allows sorting the elements to any desired order 
for editing, and facilitates greatly the process of altering a matrix for successive solutions 
under varying conditions. Restricting the transcription to one element per card also aids 
accuracy of transcription, for it leaves the transcriber free to transcribe by rows or columns 
in whatever fashion results in the least change of information from element to element. Com- 
mon information can therefore be carried from element to element on the key punch auto- 
matically. 

The most important element of accuracy control is obtained by taking row and column 
sums of sections of the matrix. These sums must be taken manually from the matrix using an 
adding machine. If the sum of the row sums agrees with the sum of the column sums, a good 
check on the accuracy of the manual summing is obtained. If the modified simplex algorithm 
is being used the coefficients of the redundant equation required to start Phase I should also be 
obtained by manual summing. After the matrix has been transcribed to punched cards, the 
cards are sorted and listed by columns and by rows, obtaining by machine the same totals that 
were taken by hand, including the grand totals. The row and column totals and the coefficients 
of the redundant equation are checked visually for correspondence. If the machine totals agree 
with the manual totals, the probability that compensating errors will occur to mask an error 
in transcription is exceedingly low. 

The work involved in taking manual totals may seem large. However, the labor 
required to insure input accuracy by other means is at least an order of magnitude greater. 
Frequently the manual totals can be taken by clerical personnel. Accounting departments have 
trained clerks who do this sort of thing routinely. 

The listing of rows and columns has another good use. Frequently it pays to examine 
at least parts of these listings to have "another look" at the structure of rows and columns 
questioning both the numerical values of the coefficients and the logic of the model. Sometimes 
it is interesting to sort out and examine meaningful groups of coefficients such as all of the 
positive and all of the negative unity coefficients, slack vectors, etc. 


MODIFICATION OF THE MATRIX 

During the development of a model, changes in the matrix are made after each trial 
solution. It is not necessary to duplicate the full editing of the matrix after each change. The 
grand sum of all the coefficients of the matrix prior to the change can be carried forward as 
a control figure. The cards for the elements deleted from this matrix should be listed and 
their total obtained. New cards being added to the matrix should similarly be listed and totaled. 
These perturbing totals should be applied to the original control sum to determine a new control 
sum. If the new matrix deck is totaled by machine and this total agrees with the new control 
figure, the chance of having an undisclosed card-handling error in the change is very low. The 
problem of making such changes in a matrix seems simple, but is not trivial by any means. 
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It is surprising how easy it is to make a card-handling error. Additional editing can beapplied 
to such an alteration of a matrix as indicated by the size and nature of the change. 4 


SUMMARY 

If the foregoing measures and others that will suggest themselves according to the 
nature of each particular model are carefully applied, transcription errors will be virtually 
eliminated and maximum progress toward an adequate model will be obtained. Considerable 
time and effort will be saved that otherwise would be required to search out error. Fewer 
trial runs will be required in the development of a model. The computing costs thus saved 
will be appreciable, but of greater importance will be the higher level of reliability that can 
be placed on the results. 
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